For this particular assignment, the data of different types of wine sales in the 20th century is to be analysed. Both of these data are from the same company but of different wines. As an analyst in the ABC Estate Wines, you are tasked to analyse and forecast Wine Sales in the 20th century.
I have generally tried to answer all the sections in the paper and have added certain parts. There are extra packages needed to fun the code which form part of my usual workflow and I will list them here as I may not be able to send an environment file. I will add the extra packages commented out in a code block to make the process simpler if you have to re-run to check anything. Because I use the tree explorer in Visual Studio I have sections and subsection in my markdowns and have tried to keep it as consistent as possible
List of extra packages and use:
Also please note I had to use cg optimisation for the SARIMA because I had an optimisation error with my LF-BFG method on my computer and my environment. I tried restarting but none of the models would converge.
#!pip install dataprep
#!pip install altair
from statsmodels.tsa.arima.model import ARIMA as ar
from statsmodels.tsa.seasonal import seasonal_decompose
from sklearn.linear_model import LinearRegression
import statsmodels as st
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
from statsmodels.tsa.arima.model import ARIMA
import itertools
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import numpy as np
import pandas as pd
from sklearn import metrics
import matplotlib.pyplot as plt
#import plotly.offline as py
%matplotlib inline
import seaborn as sns
from pylab import rcParams
from sklearn.model_selection import ParameterGrid
#setting the warnings off
import warnings
warnings.filterwarnings('ignore')
##always set the seedf for reproducibility of you machine learning working worjflow, I might not need it now but will use it
RANDOM_SEED = np.random.seed(0)
import os
###set file name to the current working directory where the full code notebook is so we can open the file
file_name = f"""{os.getcwd()}/Sparkling.csv"""
#####set the sparkling_df as and set the parse_dates to true
df_final = pd.read_csv(file_name,parse_dates=True)
#show the top values in the dataset to get a sense of the data and what it looks like
df_final.head(5)
| YearMonth | Sparkling | |
|---|---|---|
| 0 | 1980-01 | 1686 |
| 1 | 1980-02 | 1591 |
| 2 | 1980-03 | 2304 |
| 3 | 1980-04 | 1712 |
| 4 | 1980-05 | 1471 |
#show the shape of the dataset
df_final.shape
(187, 2)
#a print out of the shape of the dataset
print(f"""Your data has {df_final.shape[0]} rows and {df_final.shape[1]} columns.""")
Your data has 187 rows and 2 columns.
#information on the dataset
df_final.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 187 entries, 0 to 186 Data columns (total 2 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 YearMonth 187 non-null object 1 Sparkling 187 non-null int64 dtypes: int64(1), object(1) memory usage: 3.0+ KB
#description of the dataset statistics
df_final.describe().transpose()
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| Sparkling | 187.0 | 2402.417112 | 1295.11154 | 1070.0 | 1605.0 | 1874.0 | 2549.0 | 7242.0 |
We add a date index and do cursory checks to make sure we are running our EDA on data that makes sense and does not have:
If any of these cases are found we can treat them by interpolating for missing values and averaging if there are multiple values for the same month as the description of the data and exploration of data.head shows data at a month;y frequncy
#function that prints out whether there are missing values in the data or not
if (df_final.isnull().sum())[0] > 0 :
print(f"""There are {(df_final.isnull().sum())[0]} missing values and you have to interpolate""")
else:
print("There are no mising values in the dataset")
There are no mising values in the dataset
"""
We write code to identify duplicates in the column with the Dates, i.e. we wxpect only one record for each Year-Month value and
if there is more than one repeated value we get the average for that month
"""
duplicates = df_final[df_final["YearMonth"].duplicated()]
if not duplicates.empty:
print("There are some repeated values for some months in the column date column, need to create average:")
print(duplicates)
else:
print("No duplicates found in the date column no need to average")
No duplicates found in the date column no need to average
##creating the datetime column that is the index for the data and the start is the first date which is 1980-01-01, we set monthly frequency to the data
date_time = pd.date_range(start = '1980-01-01', periods = len(df_final) , freq = 'M')
date_time
DatetimeIndex(['1980-01-31', '1980-02-29', '1980-03-31', '1980-04-30',
'1980-05-31', '1980-06-30', '1980-07-31', '1980-08-31',
'1980-09-30', '1980-10-31',
...
'1994-10-31', '1994-11-30', '1994-12-31', '1995-01-31',
'1995-02-28', '1995-03-31', '1995-04-30', '1995-05-31',
'1995-06-30', '1995-07-31'],
dtype='datetime64[ns]', length=187, freq='M')
df_final
| YearMonth | Sparkling | |
|---|---|---|
| 0 | 1980-01 | 1686 |
| 1 | 1980-02 | 1591 |
| 2 | 1980-03 | 2304 |
| 3 | 1980-04 | 1712 |
| 4 | 1980-05 | 1471 |
| ... | ... | ... |
| 182 | 1995-03 | 1897 |
| 183 | 1995-04 | 1862 |
| 184 | 1995-05 | 1670 |
| 185 | 1995-06 | 1688 |
| 186 | 1995-07 | 2031 |
187 rows × 2 columns
#we set the date time and set the index to this column so we can set the datetime, we also drop the date_time column as it is not necessary
df_final['date_time'] = date_time
df_final.set_index(keys='date_time',inplace=True)
##removing the column YearMonth from the table
df_final.drop(['YearMonth'],axis=1,inplace=True)
df_final.head(3)
| Sparkling | |
|---|---|
| date_time | |
| 1980-01-31 | 1686 |
| 1980-02-29 | 1591 |
| 1980-03-31 | 2304 |
For the EDA I will use an automatic package to get an iverall view and then use drill downs for the the general set. My process is to look at the size of the data fist but since the data only has 187 columns I can use the package to get high level insights on the data
#this describes the statistical properties of the data
df_final.describe().transpose()
| count | mean | std | min | 25% | 50% | 75% | max | |
|---|---|---|---|---|---|---|---|---|
| Sparkling | 187.0 | 2402.417112 | 1295.11154 | 1070.0 | 1605.0 | 1874.0 | 2549.0 | 7242.0 |
The purpose of the EDA is to gain high level insights into the data and determine if the data can be used to answer the questions I have been provided with. Given that the task of the project is to gain insights on the historical sales and make accurate forecasts for future wine sales the focus of the EDA will be to do the following:
- Analysing how sales perform each year
- Analysing how sales perform each month for all the years
- The seasonal pattern of sales if any and its behavior (i.e is it additive or multiplicative in nature as well as how many seasons there are per year)
- The general trend of sales for the period provided if any is present
These insights will provide a high level understanding of the data and the direction of the analysis as well as any potential data and feature engineering techniques that need to be applied in order to run the predictions.
#install the dataprep package or module
#!pip install dataprep
from dataprep.eda import create_report
#we use the data prep create report function to create an initial EDA report on our data to give us an overall idea of what is going on
create_report(df_final.reset_index())
0%| | 0/280 [00:00<?, ?it/s]
| Number of Variables | 2 |
|---|---|
| Number of Rows | 187 |
| Missing Cells | 0 |
| Missing Cells (%) | 0.0% |
| Duplicate Rows | 0 |
| Duplicate Rows (%) | 0.0% |
| Total Size in Memory | 3.0 KB |
| Average Row Size in Memory | 16.7 B |
| Variable Types |
|
| Sparkling is skewed | Skewed |
|---|
datetime
| Distinct Count | 187.2673 |
|---|---|
| Approximate Unique (%) | 100.1% |
| Missing | 0 |
| Missing (%) | 0.0% |
| Memory Size | 1624 |
| Minimum | 1980-01-31 00:00:00 |
| Maximum | 1995-07-31 00:00:00 |
numerical
| Approximate Distinct Count | 176 |
|---|---|
| Approximate Unique (%) | 94.1% |
| Missing | 0 |
| Missing (%) | 0.0% |
| Infinite | 0 |
| Infinite (%) | 0.0% |
| Memory Size | 2992 |
| Mean | 2402.4171 |
| Minimum | 1070 |
| Maximum | 7242 |
| Zeros | 0 |
| Zeros (%) | 0.0% |
| Negatives | 0 |
| Negatives (%) | 0.0% |
| Minimum | 1070 |
|---|---|
| 5-th Percentile | 1375.6 |
| Q1 | 1605 |
| Median | 1874 |
| Q3 | 2549 |
| 95-th Percentile | 5386 |
| Maximum | 7242 |
| Range | 6172 |
| IQR | 944 |
| Mean | 2402.4171 |
|---|---|
| Standard Deviation | 1295.1115 |
| Variance | 1.6773e+06 |
| Sum | 449252 |
| Skewness | 1.803 |
| Kurtosis | 2.6055 |
| Coefficient of Variation | 0.5391 |
According to the EDA:
We will now visualise some properties of the data and sparkling series in terms of distribution of the time series
It might be helpful to plot a timeseries of our data to get a sense of the distribution of the data across time and to get a visual sense of the performance
##we visualise the time series and show the sales over time
rcParams['figure.figsize'] = 15,8
plt.plot(df_final);
plt.title("Sales of Wine")
plt.grid(False)
There looks to be clear seasonality in the data with a yearly pattern we will now further analyse the yearly trend
This gives us an understanding of the sales over the years so we can get a feel of the overall sales each year and how we are trending and performing. Since we do not have a full year for 1995, we do not include it in our Yearly EDA analysis
#add the year and month columns into the table
df_final["Year"] = df_final.index.year
df_final["month"] = df_final.index.month
df_final["Quarter"] = df_final.index.quarter
sns.relplot(
data = df_final.resample('Q').mean().reset_index(),
x = "Quarter",
y = "Sparkling",
col= "Year",
col_wrap = 4,
kind = "line"
)
<seaborn.axisgrid.FacetGrid at 0x1621d5f60>
#check that the filter excludes any data in 1995 as this is incomplete and cannot be used for a yearly analysis
df_final[:'1994'].tail(5)
| Sparkling | Year | month | Quarter | |
|---|---|---|---|---|
| date_time | ||||
| 1994-08-31 | 1495 | 1994 | 8 | 3 |
| 1994-09-30 | 2968 | 1994 | 9 | 3 |
| 1994-10-31 | 3385 | 1994 | 10 | 4 |
| 1994-11-30 | 3729 | 1994 | 11 | 4 |
| 1994-12-31 | 5999 | 1994 | 12 | 4 |
#show the performance across the years for the sparklingwine sales
# for the yearly trend we
import altair as alt
alt.Chart(df_final[:'1994'].resample('Q').mean().reset_index()).mark_line(point=True).encode(
x='quarter(date_time)',
y='Sparkling',
column='year(date_time)',
tooltip=['date_time', 'Sparkling']).properties(
title="Sales: Yearly Subseries plot",
width=100).configure_header(
titleColor='black',
titleFontSize=14,
labelColor='blue',
labelFontSize=14
)
#boxplot to see the distribution of sales in each year
rcParams['figure.figsize'] = 15,8
sns.boxplot(x = df_final[:'1994']["Year"],y = df_final[:'1994']['Sparkling'])
plt.title("Yearly Boxplot of Sales")
plt.grid();
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting. Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
#showing the yearly sales
sns.barplot(x='Year', y='Sparkling', data=df_final[:'1994'])
plt.title("Yearly Bar Chart of Sparkling Wine Sales")
plt.xlabel("Year")
plt.ylabel("Sparkling")
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting. Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Text(0, 0.5, 'Sparkling')
##showing the yearly trend, we remove 1994 as we do not have all four quarters in the data
df_final[:'1995']['Sparkling'].resample('Y').mean().plot()
plt.title("Yearly Bar Chart of Sparkling Wine Sales")
plt.xlabel("Year")
plt.ylabel("Sparkling")
Text(0, 0.5, 'Sparkling')
#looking at the yoy growth to identify any patterns in the values and only consider the data up to 1994
yearly_growth = df_final[:'1994'].groupby(df_final[:'1994'].Year)["Sparkling"].sum().pct_change()*100
yearly_growth
Year 1980 NaN 1981 -7.670915 1982 -3.454455 1983 3.392441 1984 8.598167 1985 4.252401 1986 -1.585695 1987 3.729859 1988 9.875074 1989 -5.423209 1990 -7.842763 1991 2.105118 1992 1.973840 1993 2.717842 1994 -4.540028 Name: Sparkling, dtype: float64
#plot the YOY growth or change in sales to show the performance
yearly_growth.plot(kind='line')
plt.title("YOY growth in Sparkling Wine Sales")
plt.xlabel("Year")
plt.ylabel("Sparkling")
Text(0, 0.5, 'Sparkling')
We take a look at the performance across the different quarters to understand the performance of the business within each quarter
#cplot the facet plot of each quarter
alt.Chart(df_final.resample('Q').mean().reset_index()).mark_line(point=True).encode(
x='year(date_time)' ,
y='Sparkling',
column='quarter(date_time)',
tooltip=['date_time', 'Sparkling']).properties(
title="Sales: Quarterly Subseries plot").configure_header(
titleColor='black',
titleFontSize=14,
labelColor='blue',
labelFontSize=14
)
#.properties(
# width = 'container'
#)
#boxplot of the monthly sales
sns.boxplot(x = df_final.index.month_name(),y = df_final['Sparkling'])
plt.title("Monthly Boxplot of Sales")
plt.grid();
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting. Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
#heatmap showing the distribution of the sales across quarters
sns.heatmap(pd.pivot_table(data=df_final[["Sparkling"]], index=df_final.index.year, columns=df_final.index.quarter),
square=False,
cmap='Blues',
xticklabels=["Q1", "Q2", "Q3", "Q4"]);
plt.title("Heat Map of Quarterly Sparkling Wine Sales per Year")
plt.xlabel("Quarter")
plt.ylabel("Year")
Text(158.22222222222223, 0.5, 'Year')
Here we try to ascertain the additive and multiplicative effect in each quarter and the impact each quarter has on the average sales
###get the average of the sales until the year 1994 (the one we have a full year)
avg_till_1994 = int(df_final[:'1994'].mean()[0])
##Avg sales per quarter
qrt_avg= df_final[:'1994'].groupby(df_final[:'1994'].index.quarter)["Sparkling"].mean()
##groupby quarter
qrt_table = pd.pivot_table(df_final[:'1994'][["Sparkling"]], index=df_final[:'1994'].index.quarter, columns = df_final[:'1994'].index.year)
#add qrt_avg to art_table
qrt_table["avg"] = qrt_avg
##Additive Seasonality Factor: Subtract mean from average column
qrt_table["additive"] = qrt_table["avg"] - avg_till_1994
##Multiplicative Seasonality Factor: Subtract mean from avg Column This is the amount the average for the quarter grows compared to the avergae
qrt_table["multiplicative"] = (qrt_table["avg"]/avg_till_1994).round(2)
qrt_table["mulitiplicative_percentage_change"] = (qrt_table["multiplicative"] - 1)*100
qrt_table["percentage_of_sales"] =(100* qrt_table['avg'] / qrt_table['avg'].sum()).round(2)
qrt_table.index.name = "Quarters"
qrt_table.transpose()
| Quarters | 1 | 2 | 3 | 4 | |
|---|---|---|---|---|---|
| date_time | |||||
| Sparkling | 1980 | 1860.333333 | 1520.000000 | 2134.333333 | 3954.000000 |
| 1981 | 1562.000000 | 1542.000000 | 2078.000000 | 3560.333333 | |
| 1982 | 1452.333333 | 1592.000000 | 1852.333333 | 3543.666667 | |
| 1983 | 1759.000000 | 1313.333333 | 2029.666667 | 3624.666667 | |
| 1984 | 1701.666667 | 1586.666667 | 2171.666667 | 4017.000000 | |
| 1985 | 1766.333333 | 1621.333333 | 1976.000000 | 4516.333333 | |
| 1986 | 1568.666667 | 1591.000000 | 2488.000000 | 4075.666667 | |
| 1987 | 1459.666667 | 1567.666667 | 2138.333333 | 4920.333333 | |
| 1988 | 1913.333333 | 1908.333333 | 2098.666667 | 5161.666667 | |
| 1989 | 1711.000000 | 1570.000000 | 2182.333333 | 5017.666667 | |
| 1990 | 1633.333333 | 1566.666667 | 1976.000000 | 4483.000000 | |
| 1991 | 1941.666667 | 1417.000000 | 2159.666667 | 4344.000000 | |
| 1992 | 1745.666667 | 1801.666667 | 2075.333333 | 4434.333333 | |
| 1993 | 1652.000000 | 1822.333333 | 2197.333333 | 4658.666667 | |
| 1994 | 1628.333333 | 1697.333333 | 2164.666667 | 4371.000000 | |
| avg | 1690.355556 | 1607.822222 | 2114.822222 | 4312.155556 | |
| additive | -740.644444 | -823.177778 | -316.177778 | 1881.155556 | |
| multiplicative | 0.700000 | 0.660000 | 0.870000 | 1.770000 | |
| mulitiplicative_percentage_change | -30.000000 | -34.000000 | -13.000000 | 77.000000 | |
| percentage_of_sales | 17.380000 | 16.530000 | 21.750000 | 44.340000 |
Our EDA above has already revealed that a multiplicative model may be best suited to the dataset but inorder to be certain we use both additive and multiplicative decompostion to look at the individual properties of the data
For the trend, if it is linear then the trend of the time series is additive but if it si exponentially increasing or exponential in nature then that series will be multiplicative in nature.
For the seasonality, this is caluclated relative to the level. With additive seasonality you would expect an increase to be ppositive or negative relative to the level over the period whilst when it is in terms of multiplicative seasonality you would expect that as the level changes the seasonality can vary by percentages or multiples.
##Perform the decompostion of the time series using an additive model as a test
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(df_final["Sparkling"],model='additive')
decomposition.plot();
###here we extract the individual elements of the time series and we can choose to look at the dara or not, I choose not to display it
trend = decomposition.trend
seasonality = decomposition.seasonal
residual = decomposition.resid
##looking at the residual mean
residual.mean()
-1.2088458994707845
## we plot the residuals as they are supposed to be normal accoridng to the assuomptions of the time series models
sns.distplot(residual)
<Axes: xlabel='resid', ylabel='Density'>
#we use a shapiro test to see if the additive model is actually the best fit. If it the best fit then the residuals will be normal. H0 of the Shapiro test states that the population is normally distributed
from scipy.stats import shapiro
if shapiro(residual.dropna())[1] > 0.05:
print("The distribution is normally ditributed and the time series is additive")
else:
print("The distribution is not normally ditributed and the time series is not additive")
The distribution is not normally ditributed and the time series is not additive
##Decompoisng the series into its individual components using the multiplicativbe decomposition
decomposition = seasonal_decompose(df_final["Sparkling"],model='multiplicative')
decomposition.plot();
###here we extract the individual elements of the time series and we can choose to look at the dara or not, I choose not to display it
trend = decomposition.trend
seasonality = decomposition.seasonal
residual = decomposition.resid
##looking at the residual mean
residual.mean()
0.9997456359115033
##Normality Distribution distribution
from scipy.stats import shapiro
sns.distplot(residual)
<Axes: xlabel='resid', ylabel='Density'>
#we use a shapiro test to see if the multiplicative model is actually the best fit. If it the best fit then the residuals will be normal. H0 of the Shapiro test states that the population is normally distributed
if shapiro(residual.dropna())[1] > 0.05:
print("The distribution is normally ditributed and the time series is multiplicative")
else:
print("The distribution is not normally ditributed and the time series is not multiplicative")
The distribution is normally ditributed and the time series is multiplicative
The decomposition shows the following:
We can see that a lot of the residuals for the multiplicatove series are around one and therefore we choose multiplicative decomposition
We begin our data preparation for the modeling by:
Outliers can affect how our model performs. In this case since we have identified that Q4 as higher than usual sales and that the distrbution of the data is right skewed I am expecting the outliers to be informative. We write code to identify outliers in the data and look at how we treat the outliers. From the EDA we can see that the data is right skewed showing that the sales are mostly at the lower end. We do the following to identify outliers:
##get re-accustomed with the data
df_final.head(4)
| Sparkling | Year | month | Quarter | |
|---|---|---|---|---|
| date_time | ||||
| 1980-01-31 | 1686 | 1980 | 1 | 1 |
| 1980-02-29 | 1591 | 1980 | 2 | 1 |
| 1980-03-31 | 2304 | 1980 | 3 | 1 |
| 1980-04-30 | 1712 | 1980 | 4 | 2 |
# Visual method: Box plot see if there are any outliers
sns.boxplot(x=df_final['Sparkling'])
plt.title('Box Plot of Monthly Sales')
plt.xlabel('Sparkling')
plt.show()
We can see that there is an existence of outliers (points that fall outside the "whiskers") in the data based on the boxplot on its default values. We now proceed to check for outliers using the Z-score. We have chosen the cutoff to be 3 standard deviations from the mean. So any value found over or below 3 standard deviations from the mean will be considered an outlier though we do not expect any values to be under as the distribution is right-skewed
# Statistical method: Z-score
# Calculate Z-scores
df_final['Z_score'] = np.abs((df_final['Sparkling'] - df_final['Sparkling'].mean()) /
df_final['Sparkling'].std())
# Define a threshold for outliers, usually a Z-score of 3 or -3 is considered as an outlier
#we do not need to check for the ones 3 standard deviations less than the mean since our distribution is skewed.
outliers = df_final[np.abs(df_final['Z_score']) > 3]
outliers
| Sparkling | Year | month | Quarter | Z_score | |
|---|---|---|---|---|---|
| date_time | |||||
| 1987-12-31 | 7242 | 1987 | 12 | 4 | 3.736808 |
| 1988-12-31 | 6757 | 1988 | 12 | 4 | 3.362323 |
| 1989-12-31 | 6694 | 1989 | 12 | 4 | 3.313678 |
| 1993-12-31 | 6410 | 1993 | 12 | 4 | 3.094392 |
The Z-score measures how far a data point is from the mean of a dsitribution. We have identified three 4 values with a z-score greater than 3.These points lie 3 points from the mean of the sales. The outliers are in the month of December and coincide with a large number of sales. From the inspection of the time series we would expect outliers in December which seems to imply seasonality in the sales and therefore the values contian important information that we would need to use to model.
Outliers can be treated using the following three options:
In this case we decide to clip the values so that no points lie beyond the boundary values. We caluclate our boundary values and replace the outliers with this maximum value
highest_allowed = df_final['Sparkling'].mean() + 3*df_final['Sparkling'].std()
lowest_allowed = df_final['Sparkling'].mean() - 3*df_final['Sparkling'].std()
print(f""" The highest allowed values is {highest_allowed} and the lowest allowed value is {lowest_allowed}""")
The highest allowed values is 6287.751731499824 and the lowest allowed value is -1482.9175069008938
#clipping the values in the column Sparkling to the maximum value 6287
df_final["Sparkling_Clipped"] = df_final["Sparkling"].clip(upper = 6287)
##Plot the new clipped sparkling series
rcParams['figure.figsize'] = 15,8
plt.plot(df_final["Sparkling_Clipped"]);
plt.title("Sales of Wine")
plt.ylabel('Sparkling')
plt.xlabel('Year')
plt.grid()
To train the models we split the data into training and testing datasets.
###check what df_final looks like now
df_final.head(2)
| Sparkling | Year | month | Quarter | Z_score | Sparkling_Clipped | |
|---|---|---|---|---|---|---|
| date_time | ||||||
| 1980-01-31 | 1686 | 1980 | 1 | 1 | 0.553170 | 1686 |
| 1980-02-29 | 1591 | 1980 | 2 | 1 | 0.626523 | 1591 |
#look at the unique years in the data
df_final.index.year.unique()
Int64Index([1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990,
1991, 1992, 1993, 1994, 1995],
dtype='int64', name='date_time')
We are projecting the sales for a 12 month or one year period and it is common practice for the model choices to be built on a duration longer than what we are projecting for to ensure stability of projections. In this case since 1995 is not a full year we will choose our best model based on the best performing model projecting the next 3 years and 7 months, i.e.
train : Data from 1980 until 1991 (12 years of training data) test : Data from 1990 till 1995 since (3 years 7 months)
####
# Split train-test set manually as mentioned in the notes the test data will be from 1992 till the end
train_data = df_final[df_final.index.year <= 1991]
test_data = df_final[df_final.index.year > 1991]
##we keep the original series in there for our predictions and to compare the results later
## we selected clipped and the original dataset
train_data = train_data[['Sparkling_Clipped', 'Sparkling']]
test_data = test_data[['Sparkling_Clipped', 'Sparkling']]
# Plot train-test using matplotlib.pyplot
plt.grid(False)
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("Predictions vs Actuals of Sparkling of Wine")
plt.plot(train_data['Sparkling'], 'green', label='Train data')
plt.plot(test_data['Sparkling'], 'blue', label='Test data')
plt.legend()
plt.grid();
##look at the shapes of the test and train dataset
print(train_data.shape)
print(test_data.shape)
(144, 2) (43, 2)
###look at what the train_data set looks like
train_data.head(3)
| Sparkling_Clipped | Sparkling | |
|---|---|---|
| date_time | ||
| 1980-01-31 | 1686 | 1686 |
| 1980-02-29 | 1591 | 1591 |
| 1980-03-31 | 2304 | 2304 |
We will build the following models and select which ones to use for forecasting. We build the following models:
Within each group of models we build the model and select the one with the lowest AIC.
This model will serve as our baseline model. In order to create the regression model, we create a time variable and make this the dependent variable with the sparkling sales as the independent variable.
##Get the number of rows in the data or dates in the train_data
train_data.shape[0]
144
"""
To get the time step for each individual train and test we get the shape from the time and test so that when we use predict on
the train the model knows we are on the next timestamp and the order of the time series is preserved
I prefer using the shape parameter as this allows me to potentially change my code in terms of the split percentage without having to manually find parameters
"""
train_time = [i+1 for i in range(len(train_data))]
test_time = [i+1+ train_data.shape[0] for i in range(len(test_data))]
print('Training Time instance','\n',train_time)
print('Test Time instance','\n',test_time)
Training Time instance [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144] Test Time instance [145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187]
### we create a copy of the data for the models
LinearRegression_train = train_data.copy()
LinearRegression_test = test_data.copy()
LinearRegression_train['time'] = train_time
LinearRegression_test['time'] = test_time
print('First few rows of Training Data','\n',LinearRegression_train.head(),'\n')
print('Last few rows of Training Data','\n',LinearRegression_train.tail(),'\n')
print('First few rows of Test Data','\n',LinearRegression_test.head(),'\n')
print('Last few rows of Test Data','\n',LinearRegression_test.tail(),'\n')
First few rows of Training Data
Sparkling_Clipped Sparkling time
date_time
1980-01-31 1686 1686 1
1980-02-29 1591 1591 2
1980-03-31 2304 2304 3
1980-04-30 1712 1712 4
1980-05-31 1471 1471 5
Last few rows of Training Data
Sparkling_Clipped Sparkling time
date_time
1991-08-31 1857 1857 140
1991-09-30 2408 2408 141
1991-10-31 3252 3252 142
1991-11-30 3627 3627 143
1991-12-31 6153 6153 144
First few rows of Test Data
Sparkling_Clipped Sparkling time
date_time
1992-01-31 1577 1577 145
1992-02-29 1667 1667 146
1992-03-31 1993 1993 147
1992-04-30 1997 1997 148
1992-05-31 1783 1783 149
Last few rows of Test Data
Sparkling_Clipped Sparkling time
date_time
1995-03-31 1897 1897 183
1995-04-30 1862 1862 184
1995-05-31 1670 1670 185
1995-06-30 1688 1688 186
1995-07-31 2031 2031 187
### we fit the linear model
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(LinearRegression_train[['time']],LinearRegression_train['Sparkling_Clipped'].values)
LinearRegression()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
LinearRegression()
#we run the predictions
test_predictions_model1 = lr.predict(LinearRegression_test[['time']])
LinearRegression_test['RegOnTime'] = (test_predictions_model1)
### We plot the predictions on the graph
plt.plot( train_data['Sparkling'], label='Train')
plt.plot(test_data['Sparkling'], label='Test')
plt.plot(LinearRegression_test['RegOnTime'], label='Regression On Time_Test Data')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("Linear Forecast")
plt.legend(loc='best')
plt.grid();
from sklearn import metrics
## Test Data - RMSE
rmse_model1_test = metrics.mean_squared_error(test_data['Sparkling'], (test_predictions_model1) ,squared=False)
print("For RegressionOnTime forecast on the Test Data, RMSE is %3.3f" %(rmse_model1_test))
resultsDf = pd.DataFrame({'Test RMSE': [rmse_model1_test]},index=['RegressionOnTime'])
resultsDf
For RegressionOnTime forecast on the Test Data, RMSE is 1344.645
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
The linear model performs poorly on the dataset.
In this method we assume that the value for tomorrow is the same as the one for today and so on.
##create the test and train data copies
NaiveModel_train = train_data.copy()
NaiveModel_test = test_data.copy()
#taking the last value as the "prediction" for the naive model
NaiveModel_test['naive'] = np.asarray(train_data['Sparkling_Clipped'])[len(np.asarray(train_data['Sparkling_Clipped']))-1]
NaiveModel_test['naive'].head()
date_time 1992-01-31 6153 1992-02-29 6153 1992-03-31 6153 1992-04-30 6153 1992-05-31 6153 Name: naive, dtype: int64
## we plot the data
plt.plot(NaiveModel_train['Sparkling'], label='Train')
plt.plot(test_data['Sparkling'], label='Test')
plt.plot(NaiveModel_test['naive'], label='Naive Forecast on Test Data')
plt.legend(loc='best')
plt.title("Naive Forecast")
plt.grid();
## Test Data - RMSE
rmse_model2_test = metrics.mean_squared_error(test_data['Sparkling'],NaiveModel_test['naive'],squared=False)
print("For RegressionOnTime forecast on the Test Data, RMSE is %3.3f" %(rmse_model2_test))
For RegressionOnTime forecast on the Test Data, RMSE is 3979.915
resultsDf_2 = pd.DataFrame({'Test RMSE': [rmse_model2_test]},index=['NaiveModel'])
resultsDf = pd.concat([resultsDf, resultsDf_2])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
As expected the naive model performs poorly on the data
This method uses the average of the training values
SimpleAverage_train = train_data.copy()
SimpleAverage_test = test_data.copy()
SimpleAverage_test['mean_forecast'] = train_data['Sparkling_Clipped'].mean()
SimpleAverage_test.head()
| Sparkling_Clipped | Sparkling | mean_forecast | |
|---|---|---|---|
| date_time | |||
| 1992-01-31 | 1577 | 1577 | 2396.208333 |
| 1992-02-29 | 1667 | 1667 | 2396.208333 |
| 1992-03-31 | 1993 | 1993 | 2396.208333 |
| 1992-04-30 | 1997 | 1997 | 2396.208333 |
| 1992-05-31 | 1783 | 1783 | 2396.208333 |
#plotting the predictions on one graph
plt.plot(SimpleAverage_train['Sparkling'], label='Train')
plt.plot(SimpleAverage_test['Sparkling'], label='Test')
plt.plot(SimpleAverage_test['mean_forecast'], label='Simple Average on Test Data')
plt.legend(loc='best')
plt.title("Simple Average Forecast")
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.grid();
## We Evaluate the model perfromance
rmse_model3_test = metrics.mean_squared_error(test_data['Sparkling'],SimpleAverage_test['mean_forecast'],squared=False)
print("For Simple Average forecast on the Test Data, RMSE is %3.3f" %(rmse_model3_test))
For Simple Average forecast on the Test Data, RMSE is 1268.463
resultsDf_3 = pd.DataFrame({'Test RMSE': [rmse_model3_test]},index=['SimpleAverageModel'])
resultsDf = pd.concat([resultsDf, resultsDf_3])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
Simple Exponential smoothing is a simple model as the name implies. This method assumes that a time series has the following:
According to the National Institure of Standards and technology it is often best to estimate the value ot the alpha parameter used when running the model. The alpha value determines how much prior data points influence the predction of the next value or are used to predict. A value of alpha = 1 means that prior data points enter less into the smooth.
The process will be to run the automatic fit of the model then do a grid search over variables, the automatic model will give an idea of the parameter space we can choose adn can serve as a good intuitive guide as to the model
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
##create a train and test dataset
SES_train = train_data.copy()
SES_test = test_data.copy()
model_SES = SimpleExpSmoothing(SES_train['Sparkling_Clipped'])
###fit the automated SES
model_SES_autofit = model_SES.fit(optimized=True, use_brute=True)
##print out the autofit model parameters
model_SES_autofit.params
{'smoothing_level': 0.03701387488898935,
'smoothing_trend': nan,
'smoothing_seasonal': nan,
'damping_trend': nan,
'initial_level': 1686.0,
'initial_trend': nan,
'initial_seasons': array([], dtype=float64),
'use_boxcox': False,
'lamda': None,
'remove_bias': False}
##get a summary of the model perfromance
model_SES_autofit.summary()
| Dep. Variable: | Sparkling_Clipped | No. Observations: | 144 |
|---|---|---|---|
| Model: | SimpleExpSmoothing | SSE | 234511560.409 |
| Optimized: | True | AIC | 2063.661 |
| Trend: | None | BIC | 2069.601 |
| Seasonal: | None | AICC | 2063.949 |
| Seasonal Periods: | None | Date: | Sun, 29 Oct 2023 |
| Box-Cox: | False | Time: | 18:00:33 |
| Box-Cox Coeff.: | None |
| coeff | code | optimized | |
|---|---|---|---|
| smoothing_level | 0.0370139 | alpha | True |
| initial_level | 1686.0000 | l.0 | False |
##run predicitons on the data
SES_test['auto_predict'] = model_SES_autofit.forecast(steps=len(test_data))
###take a look at our SES_test data
SES_test.head()
| Sparkling_Clipped | Sparkling | auto_predict | |
|---|---|---|---|
| date_time | |||
| 1992-01-31 | 1577 | 1577 | 2621.167381 |
| 1992-02-29 | 1667 | 1667 | 2621.167381 |
| 1992-03-31 | 1993 | 1993 | 2621.167381 |
| 1992-04-30 | 1997 | 1997 | 2621.167381 |
| 1992-05-31 | 1783 | 1783 | 2621.167381 |
## Plotting on both the Training and Test data
plt.plot(SES_train['Sparkling'], label='Train')
plt.plot(SES_test['Sparkling'], label='Test')
plt.plot((SES_test['auto_predict']), label='Alpha =0.0370 Simple Exponential Smoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("SES Auto Forecast")
plt.legend(loc='best')
plt.grid();
###Evaluating the model accuracy of the autofitted model
rmse_model_auto_SES = metrics.mean_squared_error(SES_test['Sparkling'],(SES_test['auto_predict']),squared=False)
##evaluate auto SES
resultsDf_auto_SES = pd.DataFrame({'Test RMSE': [rmse_model_auto_SES]},index=['auto_SES,Alpha =0.0370']) # Complete the code to check the perfromance of the 'Alpha = 0.05,SimpleExponentialSmoothing'
resultsDf = pd.concat([resultsDf, resultsDf_auto_SES])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
We conduct a search for the best parameters for the SES model. In this case we are searching for the alpha value
## We will now create an empty frame to pass the results in
resultsDf_SES = pd.DataFrame({'Alpha Values':[],'Train RMSE':[],'Test RMSE': []})
resultsDf_SES
| Alpha Values | Train RMSE | Test RMSE |
|---|
for i in np.arange(0.1,1,0.1):
model_SES_alpha_i = model_SES.fit(smoothing_level=i,optimized=False,use_brute=True)
SES_train['predict',i] = model_SES_alpha_i.fittedvalues
SES_test['predict',i] = model_SES_alpha_i.forecast(steps=len(SES_test))
rmse_model4_train_i = metrics.mean_squared_error(SES_train['Sparkling'],(SES_train['predict',i]),squared=False)
rmse_model4_test_i = metrics.mean_squared_error(SES_test['Sparkling'],(SES_test['predict',i]),squared=False)
resultsDf_SES = resultsDf_SES.append({'Alpha Values':i,'Train RMSE':rmse_model4_train_i
,'Test RMSE':rmse_model4_test_i}, ignore_index=True)
resultsDf_SES.sort_values(by =['Test RMSE'], ascending = True)
| Alpha Values | Train RMSE | Test RMSE | |
|---|---|---|---|
| 0 | 0.1 | 1332.559841 | 1356.185563 |
| 1 | 0.2 | 1350.607371 | 1573.611576 |
| 2 | 0.3 | 1350.051254 | 1900.040866 |
| 3 | 0.4 | 1339.105813 | 2260.068749 |
| 4 | 0.5 | 1326.469664 | 2606.296379 |
| 5 | 0.6 | 1317.231895 | 2924.118301 |
| 6 | 0.7 | 1313.280224 | 3214.744366 |
| 7 | 0.8 | 1314.955557 | 3483.731806 |
| 8 | 0.9 | 1322.275434 | 3736.981096 |
## Plotting on both the Training and Test data
plt.figure(figsize=(18,9))
plt.plot(SES_train['Sparkling'], label='Train')
plt.plot(SES_test['Sparkling'], label='Test')
plt.plot((SES_test['auto_predict']), label='Alpha =0.0370 Auto Simple Exponential Smoothing predictions on Test Set')
plt.plot((SES_test['predict', 0.1]), label='Alpha =0.1 Simple Exponential Smoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("SES Forecast")
plt.legend(loc='best')
plt.grid();
resultsDf_4 = pd.DataFrame({'Test RMSE': [resultsDf_SES.sort_values(by =['Test RMSE'], ascending = True).values[0][2]]}
,index=['Alpha=0.1,SimpleExponentialSmoothing'])
resultsDf = pd.concat([resultsDf, resultsDf_4])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
The model from the hyperparameter search performs worser than the automatic model as it has a lower RMSE. Both models however, perfrom poorly.
This model is more reliable and it does well for handling data with trends and without seasonality. The model assumes that a time series has the following:
This model accounts for alpha and beta which account for both the trend and the seasonality when modeling the data.
#creating a copy of test and training data
DES_train = train_data.copy()
DES_test = test_data.copy()
#runninf the DES AUTO MODEL
model_DES = Holt(DES_train['Sparkling_Clipped'])
##fitting the DES model
model_DES_autofit = model_DES.fit()
##Get to see the des parameters
model_DES_autofit.params
{'smoothing_level': 0.7357142857142857,
'smoothing_trend': 0.0001,
'smoothing_seasonal': nan,
'damping_trend': nan,
'initial_level': 1686.0,
'initial_trend': -95.0,
'initial_seasons': array([], dtype=float64),
'use_boxcox': False,
'lamda': None,
'remove_bias': False}
##get a summary of the model parameters
model_DES_autofit.summary()
| Dep. Variable: | Sparkling_Clipped | No. Observations: | 144 |
|---|---|---|---|
| Model: | Holt | SSE | 243084031.818 |
| Optimized: | True | AIC | 2072.831 |
| Trend: | Additive | BIC | 2084.710 |
| Seasonal: | None | AICC | 2073.444 |
| Seasonal Periods: | None | Date: | Sun, 29 Oct 2023 |
| Box-Cox: | False | Time: | 18:00:34 |
| Box-Cox Coeff.: | None |
| coeff | code | optimized | |
|---|---|---|---|
| smoothing_level | 0.7357143 | alpha | True |
| smoothing_trend | 0.0001 | beta | True |
| initial_level | 1686.0000 | l.0 | False |
| initial_trend | -95.000000 | b.0 | False |
## Prediction on the test data
DES_test['auto_predict'] = model_DES_autofit.forecast(steps=len(test_data))
DES_test.head()
| Sparkling_Clipped | Sparkling | auto_predict | |
|---|---|---|---|
| date_time | |||
| 1992-01-31 | 1577 | 1577 | 5314.317709 |
| 1992-02-29 | 1667 | 1667 | 5221.047215 |
| 1992-03-31 | 1993 | 1993 | 5127.776721 |
| 1992-04-30 | 1997 | 1997 | 5034.506227 |
| 1992-05-31 | 1783 | 1783 | 4941.235733 |
## Plotting on both the Training and Test using autofit
plt.figure(figsize=(18,9))
plt.plot(DES_train['Sparkling'], label='Train')
plt.plot(DES_test['Sparkling'], label='Test')
plt.plot((DES_test['auto_predict']), label='Alpha = 0.737,Beta = 0.0001,DoubleExponentialSmoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("DES Auto Forecast")
plt.legend(loc='best')
plt.grid();
###Evaluating the model accuracy of the autofitted model
rmse_model_auto_DES = metrics.mean_squared_error(DES_test['Sparkling'],(DES_test['auto_predict']),squared=False)
##evaluate auto model
resultsDf_auto_DES = pd.DataFrame({'Test RMSE': [rmse_model_auto_DES]},index=['auto_DES, Alpha = 0.737,Beta = 0.0001']) # Complete the code to check the perfromance of the 'Alpha = 0.05,SimpleExponentialSmoothing'
resultsDf = pd.concat([resultsDf, resultsDf_auto_DES])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
There are a few parameters that we are searching for in thsi case:
We create a parameter grid as multiple for loops would make tracing changes difficult
## we search over the parameters and use the options to try find a better model than the one given by auto predicting
parameter_grid_DES = {'smoothing_level': [.10,.20,.30,.40,.50,.60,.70,.80,.90],
'smoothing_trend':[.10,.20,.30,.40,.50,.60,.70,.80,.90],
'damping_trend': [.10,.20,.30,.40,.50,.60,.70,.80,.90],
'damped_trend': [True, False],
'exponential': [True, False] ,
#'use_boxcox':[True, False],
'remove_bias':[True, False]}
pg_DES = list(ParameterGrid(parameter_grid_DES))
##create a data frame to store the results
resultsDf_DES = pd.DataFrame(columns=['index','smoothing_level', 'smoothing_trend', 'damping_trend','damped_trend','exponential','remove_bias','Train RMSE' ,'Test RMSE' ])
####running the hyper-parameter search over values
for a,b in enumerate(pg_DES):
smoothing_level = b.get('smoothing_level')
smoothing_trend = b.get('smoothing_trend')
damping_trend = b.get('damping_trend')
damped_trend = b.get('damped_trend')
exponential = b.get('exponential')
#use_boxcox = b.get('use_boxcox')
remove_bias = b.get('remove_bias')
### we are setting new parameters so we fit a different model and not the one initialised above
model_DES_grid = Holt(DES_train['Sparkling_Clipped'], damped_trend=damping_trend,
exponential = exponential ).fit(
smoothing_level=smoothing_level,
smoothing_trend=smoothing_trend,
optimized=False,
use_brute=True,
damping_slope = damping_trend,
#use_boxcox = use_boxcox,
remove_bias = remove_bias
)
##we reference the predicitons using the enumerated number a so we can automatically reference the model
DES_train['predict', a] = model_DES_grid.fittedvalues
DES_test['predict', a] = model_DES_grid.forecast(steps=len(DES_test))
####we caluclate the error
rmse_model5_train = metrics.mean_squared_error(DES_train['Sparkling'],(DES_train['predict',a]),squared=False)
rmse_model5_test = metrics.mean_squared_error(DES_test['Sparkling'],(DES_test['predict',a]),squared=False)
#we append the results to the dataframe
resultsDf_DES = resultsDf_DES.append({'index':a,'smoothing_level':smoothing_level, 'smoothing_trend':smoothing_trend, 'damping_trend':damping_trend,
'damped_trend':damped_trend, 'exponential':exponential,'remove_bias':remove_bias, 'Train RMSE':rmse_model5_train,
'Test RMSE':rmse_model5_test}, ignore_index=True)
#we look at the sorted values in the table
resultsDf_DES.sort_values(by =['Test RMSE'], ascending = True).head()
| index | smoothing_level | smoothing_trend | damping_trend | damped_trend | exponential | remove_bias | Train RMSE | Test RMSE | |
|---|---|---|---|---|---|---|---|---|---|
| 243 | 243 | 0.1 | 0.1 | 0.1 | True | False | False | 1333.046440 | 1357.769783 |
| 3159 | 3159 | 0.1 | 0.1 | 0.1 | False | False | False | 1333.046440 | 1357.769783 |
| 81 | 81 | 0.1 | 0.1 | 0.1 | True | True | False | 1333.115094 | 1358.182750 |
| 2997 | 2997 | 0.1 | 0.1 | 0.1 | False | True | False | 1333.115094 | 1358.182750 |
| 3160 | 3160 | 0.1 | 0.2 | 0.1 | False | False | False | 1333.420423 | 1359.373198 |
Our model from the hyper-parameter search performs better than the automatic model
#we add the
resultsDf_DES_RMSE = pd.DataFrame({'Test RMSE': [resultsDf_DES.sort_values(by =['Test RMSE'], ascending = True).values[0][-1]]}
,index=['Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing']) # Complete the code to check the perfromance of the 'Alpha=0.6,Beta=0.00010,DoubleExponentialSmoothing'
resultsDf = pd.concat([resultsDf, resultsDf_DES_RMSE])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
# we get the index of the best model
index_best_model = resultsDf_DES.sort_values(by =['Test RMSE'], ascending = True).values[0][0]
print(index_best_model)
243
## Plotting on both the Training and Test data
plt.figure(figsize=(18,9))
plt.plot(DES_train['Sparkling'], label='Train')
plt.plot(DES_test['Sparkling'], label='Test')
plt.plot((DES_test['auto_predict']), label='Auto Alpha = 0.737,Beta = 0.0001,DoubleExponentialSmoothing predictions on Test Set')
plt.plot((DES_test['predict', index_best_model]), label='Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("DES Forecast")
plt.legend(loc='best')
plt.grid();
This model accounts for a time series that has the following:
Holt-Winter model is an extension of Holt's model. This takes into account all the 3 time-series components (level, trend and seasonality). The seasonality component of time series is an important component since most real-world data have a seasonal period and our own dataset has seasonality
###copy the triasning data and create a test and train setr for TES
TES_train = train_data.copy()
TES_test = test_data.copy()
###running the automated method of the TES model
model_TES = ExponentialSmoothing(TES_train['Sparkling_Clipped'],trend='additive',seasonal='multiplicative',freq='M')
#fitting the model of the TES model
model_TES_autofit = model_TES.fit()
##plotting the parameters of the TES model
model_TES_autofit.params
{'smoothing_level': 0.07701610430072522,
'smoothing_trend': 0.07701605070081585,
'smoothing_seasonal': 0.3065302280420654,
'damping_trend': nan,
'initial_level': 2352.9228257009327,
'initial_trend': -17.212966724780642,
'initial_seasons': array([0.72064637, 0.68483917, 0.89951413, 0.79977507, 0.65732154,
0.64684222, 0.87989367, 1.13119067, 0.91746882, 1.23548594,
1.90475453, 2.43819161]),
'use_boxcox': False,
'lamda': None,
'remove_bias': False}
#summary of the autofit model
model_TES_autofit.summary()
| Dep. Variable: | Sparkling_Clipped | No. Observations: | 144 |
|---|---|---|---|
| Model: | ExponentialSmoothing | SSE | 15804226.424 |
| Optimized: | True | AIC | 1703.260 |
| Trend: | Additive | BIC | 1750.777 |
| Seasonal: | Multiplicative | AICC | 1708.732 |
| Seasonal Periods: | 12 | Date: | Sun, 29 Oct 2023 |
| Box-Cox: | False | Time: | 18:01:02 |
| Box-Cox Coeff.: | None |
| coeff | code | optimized | |
|---|---|---|---|
| smoothing_level | 0.0770161 | alpha | True |
| smoothing_trend | 0.0770161 | beta | True |
| smoothing_seasonal | 0.3065302 | gamma | True |
| initial_level | 2352.9228 | l.0 | True |
| initial_trend | -17.212967 | b.0 | True |
| initial_seasons.0 | 0.7206464 | s.0 | True |
| initial_seasons.1 | 0.6848392 | s.1 | True |
| initial_seasons.2 | 0.8995141 | s.2 | True |
| initial_seasons.3 | 0.7997751 | s.3 | True |
| initial_seasons.4 | 0.6573215 | s.4 | True |
| initial_seasons.5 | 0.6468422 | s.5 | True |
| initial_seasons.6 | 0.8798937 | s.6 | True |
| initial_seasons.7 | 1.1311907 | s.7 | True |
| initial_seasons.8 | 0.9174688 | s.8 | True |
| initial_seasons.9 | 1.2354859 | s.9 | True |
| initial_seasons.10 | 1.9047545 | s.10 | True |
| initial_seasons.11 | 2.4381916 | s.11 | True |
##PREDICTION ON TEST SET
TES_test['auto_predict'] = model_TES_autofit.forecast(steps=len(test_data))
TES_test.head()
| Sparkling_Clipped | Sparkling | auto_predict | |
|---|---|---|---|
| date_time | |||
| 1992-01-31 | 1577 | 1577 | 1717.578944 |
| 1992-02-29 | 1667 | 1667 | 1606.141073 |
| 1992-03-31 | 1993 | 1993 | 1807.445563 |
| 1992-04-30 | 1997 | 1997 | 1563.969913 |
| 1992-05-31 | 1783 | 1783 | 1525.479013 |
## Plotting on both the Training and Test using autofit
plt.figure(figsize=(18,9))
plt.plot(TES_train['Sparkling'], label='Train')
plt.plot(TES_test['Sparkling'], label='Test')
plt.plot(TES_test['auto_predict'], label='Alpha=0.077016,Beta=0.077016,Gamma=0.306,Auto TripleExponentialSmoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("TES Auto Forecast")
plt.legend(loc='best')
plt.grid();
###Evaluating the model accuracy of the autofitted model
rmse_model_auto_TES = metrics.mean_squared_error(TES_test['Sparkling'],(TES_test['auto_predict']),squared=False)
##evaluate auto model
resultsDf_auto_TES = pd.DataFrame({'Test RMSE': [rmse_model_auto_TES]},index=['auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306']) # Complete the code to check the perfromance of the 'Alpha = 0.05,SimpleExponentialSmoothing'
resultsDf = pd.concat([resultsDf, resultsDf_auto_TES])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
There are a few parameters that we are searching for in thsi case:
We create a parameter grid as multiple for loops would make tracing changes difficult
##Creating the parameter grid for the TES model
param_grid_TES = {'trend': ['add', 'mul'],
'seasonal' :[ 'mul'], #we know that we have multiplicative seasonality
'seasonal_periods':[12], #we know our seasonality is 12 from the decomposition
'smoothing_level': [.10,.20,.30,.40,.50,.60,.70,.80,.90],
'smoothing_trend':[.10,.20,.30,.40,.50,.60,.70,.80,.90],
'damping_trend': [.10,.20,.30,.40,.50,.60,.70,.80,.90],
'smoothing_seasonal': [.10,.20,.30,.40,.50,.60,.70,.80,.90],
'damped_trend' : [True, False],
#'use_boxcox':[True, False],
# 'remove_bias':[True, False],
#'use_basinhopping':[True, False]
}
pg_TES = list(ParameterGrid(param_grid_TES))
##Creating the dataframe to store the results of the hyper-parameter search
df_results_TES = pd.DataFrame(columns=['index','trend','seasonal_periods','seasonal','smoothing_level', 'smoothing_trend','smoothing_seasonal',
'damped_trend','damping_trend','Train RMSE' ,'Test RMSE'])
###running the hyper-parameter search
for a,b in enumerate(pg_TES):
trend = b.get('trend')
seasonal = b.get('seasonal')
smoothing_level = b.get('smoothing_level')
seasonal_periods = b.get('seasonal_periods')
smoothing_level = b.get('smoothing_level')
smoothing_trend = b.get('smoothing_trend')
smoothing_seasonal= b.get('smoothing_seasonal')
damped_trend = b.get('damped_trend')
damping_trend = b.get('damping_trend')
#use_boxcox = b.get('use_boxcox')
#remove_bias = b.get('remove_bias')
#use_basinhopping = b.get('use_basinhopping')
#print(b)
#print(a)
##we have to initialise a different TES model from the one above as we are setting new variables in the call to the method each time
model_TES_grid = ExponentialSmoothing(TES_train['Sparkling_Clipped'], trend=trend,
damped_trend=damped_trend,
seasonal_periods=seasonal_periods, freq = 'M',
seasonal=seasonal
#,use_boxcox=use_boxcox
).fit(
# use_boxcox = use_boxcox,
smoothing_level=smoothing_level,
smoothing_trend=smoothing_trend,
smoothing_seasonal= smoothing_seasonal,
damping_trend=damping_trend,
optimized=False,
use_brute=True)
##we reference the predicitons using the enumerated number a
TES_train['predict', a] = model_TES_grid.fittedvalues
TES_test['predict', a] = model_TES_grid.forecast(steps=len(TES_test))
try:
rmse_model6_train = metrics.mean_squared_error(TES_train['Sparkling'],(TES_train['predict',a]),squared=False)
rmse_model6_test = metrics.mean_squared_error(TES_test['Sparkling'],(TES_test['predict',a]),squared=False)
df_results_TES = df_results_TES.append({'index':a, 'trend':trend, 'seasonal_periods':seasonal_periods, 'seasonal':seasonal,'smoothing_level':smoothing_level,
'smoothing_trend':smoothing_trend,'smoothing_seasonal':smoothing_seasonal, 'damped_trend':damped_trend,'damping_trend':damping_trend,
'Train RMSE':rmse_model6_train,'Test RMSE':rmse_model6_test},
ignore_index=True)
except Exception as e:
# We move on to the next one. I had this when I had a log transformation and some of the values were erroring out
#never a bad idea to have error exception handling
pass
##sort the data in the tables
df_results_TES.sort_values(by =['Test RMSE'], ascending = True).head()
#df_results_TES.to_csv('df_results_TES.csv')
| index | trend | seasonal_periods | seasonal | smoothing_level | smoothing_trend | smoothing_seasonal | damped_trend | damping_trend | Train RMSE | Test RMSE | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 17517 | 17517 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.4 | 385.434751 | 304.152583 |
| 16059 | 16059 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.3 | 385.434751 | 304.152583 |
| 13143 | 13143 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.1 | 385.434751 | 304.152583 |
| 24807 | 24807 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.9 | 385.434751 | 304.152583 |
| 18975 | 18975 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.5 | 385.434751 | 304.152583 |
##find the best model in the data
index_best_model = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][0]
print(index_best_model)
17517
The hyper-parameter search model perfroms better than the automatic TES:
The parameters for the model are:
## Plotting on both the Training and Test using autofit
plt.figure(figsize=(18,9))
plt.plot(TES_train['Sparkling'], label='Train')
plt.plot(TES_test['Sparkling'], label='Test')
#plt.plot(TES_test['auto_predict'], label='Alpha=0.0758,Beta=0.075,Gamma=0.342,TripleExponentialSmoothing predictions on Test Set')
plt.plot(TES_test['predict',index_best_model], label='Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("TES Forecast")
plt.legend(loc='best')
plt.grid();
resultsDf_DES = pd.DataFrame({'Test RMSE': [df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][-1]]}
,index=['Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing']) # Complete the code to check the perfromance of the 'Alpha=0.6,Beta=0.00010,DoubleExponentialSmoothing'
resultsDf = pd.concat([resultsDf, resultsDf_DES])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
| Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing | 304.152583 |
Given that we are trying to forecast and apply a time series model, it is important that our data be stationary, i.e. that all its statistical porpoerties (mean and the varriance) remain constant over time. Inorder for a time series to be stationary it has to have the following:
## Test for stationarity of the series - Dicky Fuller test
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(window=12).std()
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Dickey-Fuller test:
print ('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
###printing the outcome of the hypothesis test
if dftest[1] <= 0.05:
print("Reject null hypothesis and data is stationary")
else :
###
print("Fail to reject H0 and data is non-stationary ", '\n')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print (dfoutput,'\n')
test_stationarity(train_data['Sparkling_Clipped'])
Results of Dickey-Fuller Test: Fail to reject H0 and data is non-stationary Test Statistic -1.076304 p-value 0.724428 #Lags Used 13.000000 Number of Observations Used 130.000000 Critical Value (1%) -3.481682 Critical Value (5%) -2.884042 Critical Value (10%) -2.578770 dtype: float64
The time series is not stationary becasue:
##decomposing the data
decomposition = seasonal_decompose(train_data['Sparkling'],model='multiplicative')
decomposition.plot();
Since we have determined that the time series is not stationary there are a few treatments for the data. There are two major reasons for a time series being non-stationary and these tend to be:
We can choose to treat this in multiple ways:
We will start by differencing the data
"""
we test the stationarity of our data by taking the first difference and seeing if it makes a difference,
1. We create a difference time series
2. We remove-nas in place
3. We test the stationarity of the time series
"""
differenced_time_series = train_data['Sparkling_Clipped'].diff(periods=1)
differenced_time_series.dropna(inplace = True)
test_stationarity(differenced_time_series)
Results of Dickey-Fuller Test: Reject null hypothesis and data is stationary Test Statistic -7.963139e+00 p-value 2.912230e-12 #Lags Used 1.200000e+01 Number of Observations Used 1.300000e+02 Critical Value (1%) -3.481682e+00 Critical Value (5%) -2.884042e+00 Critical Value (10%) -2.578770e+00 dtype: float64
As we can see after differencing our data, the time series is now satationary because:
We can therefore proceed with forecasting and building our models in the data
we first get an idea of what values of p,d,and q would make sense for us
plot_acf(df_final["Sparkling_Clipped"].diff(periods=1).dropna())
#Plot ACF and PACF using statsmodels
plot_pacf(df_final["Sparkling_Clipped"].diff(periods=1).dropna());
Checking the order of differencing we can use
result_1d = adfuller(df_final["Sparkling"].diff().dropna())
print(f""" p-value {result_1d[1]}""")
result_2d= adfuller(df_final["Sparkling"].diff(periods=2).dropna())
print(f""" p-value {result_2d[1]}""")
p-value 0.0 p-value 0.0
In both plots we can see that the first lag is the most significant and we therefore consider p to be 1. To determine q we can look at the ACF plot and consider at the number of lags crossing the threshold as this helps us determine how much of the past data would be significant enough to consider for the future. In this case we would want the ones with high correlation as they would contribute the most and are the most useful in predicting our future values. We would expect our q to be 2 to 3. We can see from the AD-fuller test above that both the first and second order differencing produces a stationary time series so we expect a value of 1 for d to be sufficient. This analysis will help us set a range for the parameter searches.
Our values for the parameters for grid searches are therefore going to be:
###Creating the parameters for the models
import itertools
p = q = range(0,4)
d = range(1,2)
pdq = list(itertools.product(p,d,q))
print('Examples of the parameter combinations for the models')
for i in range(0,len(pdq)):
print('Model : {}'.format(pdq[i]))
Examples of the parameter combinations for the models Model : (0, 1, 0) Model : (0, 1, 1) Model : (0, 1, 2) Model : (0, 1, 3) Model : (1, 1, 0) Model : (1, 1, 1) Model : (1, 1, 2) Model : (1, 1, 3) Model : (2, 1, 0) Model : (2, 1, 1) Model : (2, 1, 2) Model : (2, 1, 3) Model : (3, 1, 0) Model : (3, 1, 1) Model : (3, 1, 2) Model : (3, 1, 3)
# Creating an empty Dataframe with column names only
ARIMA_AIC = pd.DataFrame(columns=['param', 'AIC'])
ARIMA_AIC
| param | AIC |
|---|
###run hyper - parameter search
from statsmodels.tsa.arima.model import ARIMA
for param in pdq: # running a loop within the pdq parameters defined by itertools
ARIMA_model = ARIMA(train_data['Sparkling_Clipped'].values, order=param).fit()
print('ARIMA{} - AIC{}'.format(param, ARIMA_model.aic))
#printing the parameters and the AIC from the fitted models
ARIMA_AIC = ARIMA_AIC.append({'param': param, 'AIC': ARIMA_model.aic},ignore_index=True)
#appending the AIC values and the model parameters to the previously created data frame
#for easier understanding and sorting of the AIC values
ARIMA(0, 1, 0) - AIC2462.923400148405 ARIMA(0, 1, 1) - AIC2459.1483109544283 ARIMA(0, 1, 2) - AIC2427.157170600601 ARIMA(0, 1, 3) - AIC2426.2964726118053 ARIMA(1, 1, 0) - AIC2462.176371753829 ARIMA(1, 1, 1) - AIC2428.574792704891 ARIMA(1, 1, 2) - AIC2427.061431733311 ARIMA(1, 1, 3) - AIC2428.027348160358 ARIMA(2, 1, 0) - AIC2455.6158026664075 ARIMA(2, 1, 1) - AIC2426.0100062824467 ARIMA(2, 1, 2) - AIC2404.1798175402737 ARIMA(2, 1, 3) - AIC2424.3235635933197 ARIMA(3, 1, 0) - AIC2453.223330622951 ARIMA(3, 1, 1) - AIC2427.62908683188 ARIMA(3, 1, 2) - AIC2422.586082662382 ARIMA(3, 1, 3) - AIC2412.311062543722
##see the top perfroming model
ARIMA_AIC.sort_values(by='AIC',ascending=True).head()
| param | AIC | |
|---|---|---|
| 10 | (2, 1, 2) | 2404.179818 |
| 15 | (3, 1, 3) | 2412.311063 |
| 14 | (3, 1, 2) | 2422.586083 |
| 11 | (2, 1, 3) | 2424.323564 |
| 9 | (2, 1, 1) | 2426.010006 |
###choose the model that has performed the best as use the parameters
auto_ARIMA = ARIMA(train_data['Sparkling_Clipped'], order=(2,1,2))
results_auto_ARIMA = auto_ARIMA.fit()
print(results_auto_ARIMA.summary())
SARIMAX Results
==============================================================================
Dep. Variable: Sparkling_Clipped No. Observations: 144
Model: ARIMA(2, 1, 2) Log Likelihood -1197.090
Date: Sun, 29 Oct 2023 AIC 2404.180
Time: 18:05:20 BIC 2418.994
Sample: 01-31-1980 HQIC 2410.200
- 12-31-1991
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 1.3289 0.043 30.646 0.000 1.244 1.414
ar.L2 -0.5729 0.063 -9.122 0.000 -0.696 -0.450
ma.L1 -1.9919 0.112 -17.746 0.000 -2.212 -1.772
ma.L2 0.9999 0.113 8.875 0.000 0.779 1.221
sigma2 1.021e+06 2.21e-07 4.62e+12 0.000 1.02e+06 1.02e+06
===================================================================================
Ljung-Box (L1) (Q): 0.08 Jarque-Bera (JB): 7.32
Prob(Q): 0.77 Prob(JB): 0.03
Heteroskedasticity (H): 1.72 Skew: 0.49
Prob(H) (two-sided): 0.06 Kurtosis: 3.53
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 7.03e+27. Standard errors may be unstable.
results_auto_ARIMA.plot_diagnostics();
Most of the plots are as expected but the QQ plot shows a lot of outliers at the tail end of the model
predicted_auto_ARIMA = results_auto_ARIMA.forecast(steps=len(test_data))
## Mean Absolute Percentage Error (MAPE) - Function Definition
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean((np.abs(y_true-y_pred))/(y_true))*100
## Importing the mean_squared_error function from sklearn to calculate the RMSE
from sklearn.metrics import mean_squared_error
rmse = mean_squared_error(test_data['Sparkling'],(predicted_auto_ARIMA),squared=False)
mape = mean_absolute_percentage_error(test_data['Sparkling'],(predicted_auto_ARIMA))
print('RMSE:',rmse,'\nMAPE:',mape)
RMSE: 1300.8578710583665 MAPE: 43.382776272537896
from math import sqrt
from sklearn.metrics import mean_squared_error
rmse = sqrt(mean_squared_error(test_data['Sparkling'],(predicted_auto_ARIMA)))
print(rmse)
1300.8578710583665
ARIMA_AIC.sort_values(by='AIC',ascending=True)
| param | AIC | |
|---|---|---|
| 10 | (2, 1, 2) | 2404.179818 |
| 15 | (3, 1, 3) | 2412.311063 |
| 14 | (3, 1, 2) | 2422.586083 |
| 11 | (2, 1, 3) | 2424.323564 |
| 9 | (2, 1, 1) | 2426.010006 |
| 3 | (0, 1, 3) | 2426.296473 |
| 6 | (1, 1, 2) | 2427.061432 |
| 2 | (0, 1, 2) | 2427.157171 |
| 13 | (3, 1, 1) | 2427.629087 |
| 7 | (1, 1, 3) | 2428.027348 |
| 5 | (1, 1, 1) | 2428.574793 |
| 12 | (3, 1, 0) | 2453.223331 |
| 8 | (2, 1, 0) | 2455.615803 |
| 1 | (0, 1, 1) | 2459.148311 |
| 4 | (1, 1, 0) | 2462.176372 |
| 0 | (0, 1, 0) | 2462.923400 |
resultsDf_ARIMA = pd.DataFrame({'Test RMSE': [sqrt(mean_squared_error(test_data['Sparkling'],(predicted_auto_ARIMA)))]}
,index=['ARIMA(2,1,2)'])
resultsDf = pd.concat([resultsDf, resultsDf_ARIMA])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
| Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing | 304.152583 |
| ARIMA(2,1,2) | 1300.857871 |
resultsDf_ARIMA_singular = pd.DataFrame({'RMSE': rmse,'MAPE':mape}
,index=['ARIMA(2,1,2)'])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
| Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing | 304.152583 |
| ARIMA(2,1,2) | 1300.857871 |
plt.plot(train_data['Sparkling'], label='Train')
plt.plot(test_data['Sparkling'], label='Test')
plt.plot((predicted_auto_ARIMA), label='ARIMA(2,1,2) predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("ARIMA Forecast")
plt.legend(loc='best')
plt.grid();
We begin by getting an understanding of our seasonality over time. I will choose the time horizon of the test set and see what the seasonality is over that period.
#Plot ACF and PACF using statsmodels
plot_acf(df_final["Sparkling"].diff(periods=1).dropna(), lags = len(test_data));
We can see a seasonality of 12 in the ACF plot, the decompostion in our EDA shows that the seasonal period is also 12 so we will use 12 as our seasonal period.
##Creating the model combinations
import itertools
p = q = range(0, 4)
d= range(1,2)
D = range(0,1)
pdq = list(itertools.product(p, d, q))
PDQ = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, D, q))]
print('Examples of the parameter combinations for the Model are')
for i in range(1,len(pdq)):
print('Model: {}{}'.format(pdq[i], PDQ[i]))
Examples of the parameter combinations for the Model are Model: (0, 1, 1)(0, 0, 1, 12) Model: (0, 1, 2)(0, 0, 2, 12) Model: (0, 1, 3)(0, 0, 3, 12) Model: (1, 1, 0)(1, 0, 0, 12) Model: (1, 1, 1)(1, 0, 1, 12) Model: (1, 1, 2)(1, 0, 2, 12) Model: (1, 1, 3)(1, 0, 3, 12) Model: (2, 1, 0)(2, 0, 0, 12) Model: (2, 1, 1)(2, 0, 1, 12) Model: (2, 1, 2)(2, 0, 2, 12) Model: (2, 1, 3)(2, 0, 3, 12) Model: (3, 1, 0)(3, 0, 0, 12) Model: (3, 1, 1)(3, 0, 1, 12) Model: (3, 1, 2)(3, 0, 2, 12) Model: (3, 1, 3)(3, 0, 3, 12)
##creating the empty dataframe to plot the results in
SARIMA_AIC = pd.DataFrame(columns=['param','seasonal', 'AIC'])
SARIMA_AIC
| param | seasonal | AIC |
|---|
train_data.head()
| Sparkling_Clipped | Sparkling | |
|---|---|---|
| date_time | ||
| 1980-01-31 | 1686 | 1686 |
| 1980-02-29 | 1591 | 1591 |
| 1980-03-31 | 2304 | 2304 |
| 1980-04-30 | 1712 | 1712 |
| 1980-05-31 | 1471 | 1471 |
#using threading because my lfgbs optimizer stopped working a few hours before submission
#please note this part takes time
import warnings
import statsmodels.api as sm
import pandas as pd
import concurrent.futures
warnings.filterwarnings('ignore')
# Assuming pdq and PDQ are lists of parameters and train_data is a DataFrame
SARIMA_AIC = pd.DataFrame(columns=['param', 'seasonal', 'AIC'])
def fit_sarima_model(param, param_seasonal):
SARIMA_model = sm.tsa.statespace.SARIMAX(train_data['Sparkling_Clipped'].values,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results_SARIMA = SARIMA_model.fit(maxiter=200, method='cg')
#results_SARIMA = SARIMA_model.fit(maxiter=800)
print(f'SARIMA{param}x{param_seasonal} - AIC:{results_SARIMA.aic}')
return {'param': param, 'seasonal': param_seasonal, 'AIC': results_SARIMA.aic}
# Using ThreadPoolExecutor to parallelize the execution
with concurrent.futures.ThreadPoolExecutor() as executor:
futures = [executor.submit(fit_sarima_model, param, param_seasonal)
for param in pdq for param_seasonal in PDQ]
for future in concurrent.futures.as_completed(futures):
SARIMA_AIC = SARIMA_AIC.append(future.result(), ignore_index=True)
Optimization terminated successfully.
Current function value: 7.751307
Iterations: 7
Function evaluations: 24
Gradient evaluations: 24
Optimization terminated successfully.
Current function value: 9.682606
Iterations: 6
Function evaluations: 73
Gradient evaluations: 73
SARIMA(0, 1, 0)x(0, 0, 0, 12) - AIC:3887.1436085126074
SARIMA(0, 1, 0)x(0, 0, 1, 12) - AIC:3041.954643080705
Current function value: 25.913435
Iterations: 0
Function evaluations: 18
Gradient evaluations: 6
Optimization terminated successfully.
Current function value: 6.729041
Iterations: 9
Function evaluations: 35
Gradient evaluations: 35
SARIMA(0, 1, 0)x(0, 0, 3, 12) - AIC:8235.438067536941
SARIMA(0, 1, 0)x(0, 0, 2, 12) - AIC:2617.069292460136
Optimization terminated successfully.
Current function value: 6.958860
Iterations: 29
Function evaluations: 72
Gradient evaluations: 72
SARIMA(0, 1, 0)x(1, 0, 0, 12) - AIC:3509.674459141471
Optimization terminated successfully.
Current function value: 8.427763
Iterations: 6
Function evaluations: 22
Gradient evaluations: 22
SARIMA(0, 1, 1)x(0, 0, 0, 12) - AIC:3301.721371587172
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 12
Gradient evaluations: 12
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 11
Gradient evaluations: 11
SARIMA(0, 1, 0)x(2, 0, 3, 12) - AIC:7487.259868847646SARIMA(0, 1, 0)x(1, 0, 3, 12) - AIC:8093.113009771254
Optimization terminated successfully.
Current function value: 6.340852
Iterations: 25
Function evaluations: 70
Gradient evaluations: 70
SARIMA(0, 1, 0)x(2, 0, 0, 12) - AIC:3062.545945612424
Optimization terminated successfully.
Current function value: 6.939680
Iterations: 16
Function evaluations: 38
Gradient evaluations: 36
SARIMA(0, 1, 1)x(1, 0, 0, 12) - AIC:2811.119296644377
Optimization terminated successfully.
Current function value: 6.890972
Iterations: 41
Function evaluations: 102
Gradient evaluations: 102
Optimization terminated successfully.
Current function value: 7.358619
Iterations: 20
Function evaluations: 53
Gradient evaluations: 53
SARIMA(0, 1, 0)x(1, 0, 1, 12) - AIC:2804.33798370544
SARIMA(0, 1, 1)x(0, 0, 1, 12) - AIC:2863.830643987012
Optimization terminated successfully.
Current function value: 6.287120
Iterations: 35
Function evaluations: 76
Gradient evaluations: 76
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 11
Gradient evaluations: 11
SARIMA(0, 1, 0)x(2, 0, 1, 12) - AIC:2440.413640333953
SARIMA(0, 1, 0)x(3, 0, 3, 12) - AIC:7566.197506320239
Optimization terminated successfully.
Current function value: 6.229900
Iterations: 31
Function evaluations: 84
Gradient evaluations: 84
SARIMA(0, 1, 0)x(1, 0, 2, 12) - AIC:2420.8491108170383
Current function value: 23.493788
Iterations: 1
Function evaluations: 37
Gradient evaluations: 24
SARIMA(0, 1, 1)x(0, 0, 3, 12) - AIC:8158.871778693826
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 1)x(1, 0, 3, 12) - AIC:7727.12144020455
Optimization terminated successfully.
Current function value: 6.270336
Iterations: 27
Function evaluations: 60
Gradient evaluations: 59
SARIMA(0, 1, 1)x(2, 0, 0, 12) - AIC:2546.199154769182
Optimization terminated successfully.
Current function value: 6.640792
Iterations: 41
Function evaluations: 110
Gradient evaluations: 109
SARIMA(0, 1, 1)x(1, 0, 1, 12) - AIC:2710.148681449255
Optimization terminated successfully.
Current function value: 6.531468
Iterations: 39
Function evaluations: 116
Gradient evaluations: 116
SARIMA(0, 1, 1)x(0, 0, 2, 12) - AIC:2511.5669221735793
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 1)x(2, 0, 3, 12) - AIC:7760.205212457011
Optimization terminated successfully.
Current function value: 8.356468
Iterations: 11
Function evaluations: 49
Gradient evaluations: 49
SARIMA(0, 1, 2)x(0, 0, 0, 12) - AIC:3275.5358529353794
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 1)x(3, 0, 3, 12) - AIC:7824.555924176051
Optimization terminated successfully.
Current function value: 7.290458
Iterations: 29
Function evaluations: 84
Gradient evaluations: 84
SARIMA(0, 1, 2)x(0, 0, 1, 12) - AIC:2840.5994191671143
Optimization terminated successfully.
Current function value: 5.640160
Iterations: 28
Function evaluations: 62
Gradient evaluations: 62
SARIMA(0, 1, 1)x(3, 0, 0, 12) - AIC:2292.7920280731914
Optimization terminated successfully.
Current function value: 6.803916
Iterations: 62
Function evaluations: 151
Gradient evaluations: 151
SARIMA(0, 1, 2)x(1, 0, 0, 12) - AIC:2767.911131391426
Optimization terminated successfully.
Current function value: 5.662339
Iterations: 100
Function evaluations: 254
Gradient evaluations: 254
SARIMA(0, 1, 0)x(3, 0, 0, 12) - AIC:2890.6371489075823
Current function value: 23.240226
Iterations: 1
Function evaluations: 52
Gradient evaluations: 40
SARIMA(0, 1, 2)x(0, 0, 3, 12) - AIC:8085.042354818635
Optimization terminated successfully.
Current function value: 6.229713
Iterations: 142
Function evaluations: 338
Gradient evaluations: 338
SARIMA(0, 1, 0)x(2, 0, 2, 12) - AIC:2418.6042057746436
Optimization terminated successfully.
Current function value: 5.987934
Iterations: 108
Function evaluations: 245
Gradient evaluations: 245
SARIMA(0, 1, 1)x(1, 0, 2, 12) - AIC:2453.5194097863673
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 2)x(1, 0, 3, 12) - AIC:7679.178520952888
Optimization terminated successfully.
Current function value: 6.469502
Iterations: 71
Function evaluations: 169
Gradient evaluations: 169
SARIMA(0, 1, 2)x(0, 0, 2, 12) - AIC:2490.265661257791
Optimization terminated successfully.
Current function value: 6.102234
Iterations: 172
Function evaluations: 377
Gradient evaluations: 377
SARIMA(0, 1, 1)x(2, 0, 1, 12) - AIC:2499.4572061193203
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 2)x(2, 0, 3, 12) - AIC:7711.966257279268
Optimization terminated successfully.
Current function value: 6.127441
Iterations: 76
Function evaluations: 158
Gradient evaluations: 158
SARIMA(0, 1, 2)x(2, 0, 0, 12) - AIC:2508.450955395183
Optimization terminated successfully.
Current function value: 5.657563
Iterations: 102
Function evaluations: 218
Gradient evaluations: 218
SARIMA(0, 1, 0)x(3, 0, 1, 12) - AIC:2196.0806962769116
Optimization terminated successfully.
Current function value: 6.535281
Iterations: 182
Function evaluations: 409
Gradient evaluations: 409
SARIMA(0, 1, 2)x(1, 0, 1, 12) - AIC:2648.6018710796207
Optimization terminated successfully.
Current function value: 8.179172
Iterations: 12
Function evaluations: 65
Gradient evaluations: 65
SARIMA(0, 1, 3)x(0, 0, 0, 12) - AIC:3210.6746550658863Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 2)x(3, 0, 3, 12) - AIC:7775.7231526960895
Optimization terminated successfully.
Current function value: 5.987897
Iterations: 190
Function evaluations: 408
Gradient evaluations: 408
SARIMA(0, 1, 1)x(2, 0, 2, 12) - AIC:2455.31254999365
Current function value: 22.734271
Iterations: 1
Function evaluations: 21
Gradient evaluations: 9
SARIMA(0, 1, 3)x(0, 0, 3, 12) - AIC:7221.965069751158
Optimization terminated successfully.
Current function value: 7.113263
Iterations: 71
Function evaluations: 177
Gradient evaluations: 177
SARIMA(0, 1, 3)x(0, 0, 1, 12) - AIC:2778.5170915574126
Optimization terminated successfully.
Current function value: 6.747850
Iterations: 52
Function evaluations: 122
Gradient evaluations: 122
SARIMA(0, 1, 3)x(1, 0, 0, 12) - AIC:2746.5996227307014
Optimization terminated successfully.
Current function value: 6.292769
Iterations: 50
Function evaluations: 121
Gradient evaluations: 121
SARIMA(0, 1, 3)x(0, 0, 2, 12) - AIC:2431.8940607574064
Current function value: 6.204830
Iterations: 200
Function evaluations: 412
Gradient evaluations: 412
Optimization terminated successfully.
Current function value: 5.502948
Iterations: 69
Function evaluations: 149
Gradient evaluations: 149
SARIMA(0, 1, 2)x(2, 0, 1, 12) - AIC:2520.773774214859
SARIMA(0, 1, 2)x(3, 0, 0, 12) - AIC:2256.9300304121434
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 3)x(1, 0, 3, 12) - AIC:7689.741293759455
Optimization terminated successfully.
Current function value: 5.657521
Iterations: 138
Function evaluations: 311
Gradient evaluations: 311
SARIMA(0, 1, 0)x(3, 0, 2, 12) - AIC:2199.1387840798607
Current function value: 5.951340
Iterations: 200
Function evaluations: 439
Gradient evaluations: 439
SARIMA(0, 1, 2)x(1, 0, 2, 12) - AIC:2410.2479276862637
Optimization terminated successfully.
Current function value: 6.453198
Iterations: 134
Function evaluations: 300
Gradient evaluations: 300
SARIMA(0, 1, 3)x(1, 0, 1, 12) - AIC:2618.3113442755957
Optimization terminated successfully.
Current function value: 5.487697
Iterations: 138
Function evaluations: 294
Gradient evaluations: 294
SARIMA(0, 1, 1)x(3, 0, 1, 12) - AIC:2251.7899534305884
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 3)x(2, 0, 3, 12) - AIC:7722.2329941600865
Current function value: 5.944411
Iterations: 200
Function evaluations: 423
Gradient evaluations: 423
SARIMA(0, 1, 2)x(2, 0, 2, 12) - AIC:2413.8486259725596
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(0, 1, 3)x(3, 0, 3, 12) - AIC:7785.396073274355
Optimization terminated successfully.
Current function value: 6.088331
Iterations: 96
Function evaluations: 204
Gradient evaluations: 204
SARIMA(0, 1, 3)x(2, 0, 0, 12) - AIC:2489.9634405186316
Optimization terminated successfully.
Current function value: 8.554006
Iterations: 23
Function evaluations: 86
Gradient evaluations: 86
SARIMA(1, 1, 0)x(0, 0, 0, 12) - AIC:4199.538753322784
Optimization terminated successfully.
Current function value: 7.739127
Iterations: 23
Function evaluations: 61
Gradient evaluations: 61
SARIMA(1, 1, 0)x(0, 0, 1, 12) - AIC:3065.769619362932
Optimization terminated successfully.
Current function value: 6.719620
Iterations: 24
Function evaluations: 63
Gradient evaluations: 63
SARIMA(1, 1, 0)x(0, 0, 2, 12) - AIC:2610.3627191844475
Optimization terminated successfully.
Current function value: 6.865454
Iterations: 87
Function evaluations: 200
Gradient evaluations: 200
SARIMA(1, 1, 0)x(1, 0, 0, 12) - AIC:3447.4866976155354
Optimization terminated successfully.
Current function value: 6.843516
Iterations: 47
Function evaluations: 102
Gradient evaluations: 102
SARIMA(1, 1, 0)x(1, 0, 1, 12) - AIC:2779.519946045093
Current function value: 25.925518
Iterations: 1
Function evaluations: 61
Gradient evaluations: 48
SARIMA(1, 1, 0)x(0, 0, 3, 12) - AIC:8240.146956186698
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 0)x(1, 0, 3, 12) - AIC:7805.315908812001
Optimization terminated successfully.
Current function value: 6.071349
Iterations: 135
Function evaluations: 306
Gradient evaluations: 306
SARIMA(0, 1, 3)x(2, 0, 1, 12) - AIC:2486.2526342006176
Optimization terminated successfully.
Current function value: 6.175653
Iterations: 54
Function evaluations: 124
Gradient evaluations: 124
SARIMA(1, 1, 0)x(1, 0, 2, 12) - AIC:2515.5516100217737
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 6.186875
Iterations: 51
Function evaluations: 109
Gradient evaluations: 109
SARIMA(1, 1, 0)x(2, 0, 3, 12) - AIC:7838.695716990902
SARIMA(1, 1, 0)x(2, 0, 1, 12) - AIC:2520.453539994749
Optimization terminated successfully.
Current function value: 6.507488
Iterations: 142
Function evaluations: 318
Gradient evaluations: 318
SARIMA(1, 1, 0)x(2, 0, 0, 12) - AIC:2879.111185227929
Current function value: 5.486936
Iterations: 200
Function evaluations: 399
Gradient evaluations: 399
SARIMA(0, 1, 1)x(3, 0, 2, 12) - AIC:2253.5693203586366
Optimization terminated successfully.
Current function value: 5.497469
Iterations: 149
Function evaluations: 319
Gradient evaluations: 319
SARIMA(0, 1, 2)x(3, 0, 1, 12) - AIC:2256.376338599889
Optimization terminated successfully.
Current function value: 5.474050
Iterations: 87
Function evaluations: 192
Gradient evaluations: 192
Current function value: 5.870026
Iterations: 200
Function evaluations: 422
Gradient evaluations: 422
SARIMA(0, 1, 3)x(3, 0, 0, 12) - AIC:2242.171316304141
Optimization terminated successfully.
Current function value: 8.420898
Iterations: 21
Function evaluations: 46
Gradient evaluations: 46
SARIMA(1, 1, 1)x(0, 0, 0, 12) - AIC:3298.6543715474377
SARIMA(0, 1, 3)x(1, 0, 2, 12) - AIC:2381.038366036927
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 0)x(3, 0, 3, 12) - AIC:7903.640245011804
Current function value: 23.488527
Iterations: 0
Function evaluations: 7
Gradient evaluations: 5
SARIMA(1, 1, 1)x(0, 0, 3, 12) - AIC:8167.707946795989
Optimization terminated successfully.
Current function value: 7.352753
Iterations: 38
Function evaluations: 79
Gradient evaluations: 79
SARIMA(1, 1, 1)x(0, 0, 1, 12) - AIC:2861.960722918811
Optimization terminated successfully.
Current function value: 6.858622
Iterations: 43
Function evaluations: 76
Gradient evaluations: 76
SARIMA(1, 1, 1)x(1, 0, 0, 12) - AIC:2782.127228087847
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 6.526047
Iterations: 40
Function evaluations: 90
Gradient evaluations: 90
SARIMA(1, 1, 1)x(1, 0, 3, 12) - AIC:7780.627003442225
SARIMA(1, 1, 1)x(0, 0, 2, 12) - AIC:2510.098880777272
Optimization terminated successfully.
Current function value: 6.175112
Iterations: 163
Function evaluations: 368
Gradient evaluations: 368
SARIMA(1, 1, 0)x(2, 0, 2, 12) - AIC:2518.1601912365854
Optimization terminated successfully.
Current function value: 6.706098
Iterations: 77
Function evaluations: 158
Gradient evaluations: 158
SARIMA(1, 1, 1)x(1, 0, 1, 12) - AIC:2730.8149325172494
Optimization terminated successfully.
Current function value: 6.185196
Iterations: 40
Function evaluations: 78
Gradient evaluations: 78
SARIMA(1, 1, 1)x(2, 0, 0, 12) - AIC:2516.2823227819595
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 1)x(2, 0, 3, 12) - AIC:7813.710775692123
Current function value: 5.863493
Iterations: 200
Function evaluations: 475
Gradient evaluations: 475
SARIMA(0, 1, 3)x(2, 0, 2, 12) - AIC:2385.50290582914
Optimization terminated successfully.
Current function value: 5.556197
Iterations: 34
Function evaluations: 71
Gradient evaluations: 71
SARIMA(1, 1, 1)x(3, 0, 0, 12) - AIC:2263.314681550093
Optimization terminated successfully.
Current function value: 5.562991
Iterations: 95
Function evaluations: 205
Gradient evaluations: 205
SARIMA(1, 1, 0)x(3, 0, 1, 12) - AIC:2269.2661096705747
Current function value: 9.707374
Iterations: 200
Function evaluations: 406
Gradient evaluations: 406
SARIMA(1, 1, 0)x(3, 0, 0, 12) - AIC:2234.179850107131
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 8.243157
Iterations: 49
Function evaluations: 143
Gradient evaluations: 143
SARIMA(1, 1, 2)x(0, 0, 0, 12) - AIC:3241.498435691376
SARIMA(1, 1, 1)x(3, 0, 3, 12) - AIC:7878.061487413748
Current function value: 23.485586
Iterations: 1
Function evaluations: 19
Gradient evaluations: 10
SARIMA(1, 1, 2)x(0, 0, 3, 12) - AIC:7463.465966731997
Optimization terminated successfully.
Current function value: 7.187984
Iterations: 64
Function evaluations: 149
Gradient evaluations: 149
SARIMA(1, 1, 2)x(0, 0, 1, 12) - AIC:2811.7400660805843
Optimization terminated successfully.
Current function value: 5.980976
Iterations: 136
Function evaluations: 299
Gradient evaluations: 299
SARIMA(1, 1, 1)x(1, 0, 2, 12) - AIC:2451.1531023031857
Optimization terminated successfully.
Current function value: 6.705767
Iterations: 96
Function evaluations: 183
Gradient evaluations: 183
SARIMA(1, 1, 2)x(1, 0, 0, 12) - AIC:2742.001849580438
Optimization terminated successfully.
Current function value: 6.041225
Iterations: 149
Function evaluations: 338
Gradient evaluations: 338
Current function value: 5.503325
Iterations: 200
Function evaluations: 410
Gradient evaluations: 410
SARIMA(1, 1, 1)x(2, 0, 1, 12) - AIC:2475.224466441966
SARIMA(0, 1, 2)x(3, 0, 2, 12) - AIC:2232.4565789361563
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 2)x(1, 0, 3, 12) - AIC:7682.402027821568
Optimization terminated successfully.
Current function value: 6.363932
Iterations: 110
Function evaluations: 230
Gradient evaluations: 230
SARIMA(1, 1, 2)x(0, 0, 2, 12) - AIC:2548.846951068069
Optimization terminated successfully.
Current function value: 6.512554
Iterations: 145
Function evaluations: 313
Gradient evaluations: 313
SARIMA(1, 1, 2)x(1, 0, 1, 12) - AIC:2643.828477037553
Optimization terminated successfully.
Current function value: 6.100146
Iterations: 164
Function evaluations: 335
Gradient evaluations: 335
SARIMA(1, 1, 1)x(2, 0, 2, 12) - AIC:2471.5628162499947
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 2)x(2, 0, 3, 12) - AIC:7715.189764149396
Optimization terminated successfully.
Current function value: 6.049810
Iterations: 93
Function evaluations: 198
Gradient evaluations: 198
SARIMA(1, 1, 2)x(2, 0, 0, 12) - AIC:2480.03661359743
Optimization terminated successfully.
Current function value: 5.562990
Iterations: 139
Function evaluations: 307
Gradient evaluations: 307
SARIMA(1, 1, 0)x(3, 0, 2, 12) - AIC:2271.2543802461214
Optimization terminated successfully.
Current function value: 5.555745
Iterations: 86
Function evaluations: 182
Gradient evaluations: 182
SARIMA(1, 1, 1)x(3, 0, 1, 12) - AIC:2265.1359502337605
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 2)x(3, 0, 3, 12) - AIC:7778.946659564769
Optimization terminated successfully.
Current function value: 5.472738
Iterations: 179
Function evaluations: 389
Gradient evaluations: 389
SARIMA(0, 1, 3)x(3, 0, 1, 12) - AIC:2243.601382810842
Current function value: 8.172270
Iterations: 200
Function evaluations: 431
Gradient evaluations: 431
SARIMA(1, 1, 3)x(0, 0, 0, 12) - AIC:3209.956590610387
Optimization terminated successfully.
Current function value: 7.109882
Iterations: 127
Function evaluations: 268
Gradient evaluations: 268
SARIMA(1, 1, 3)x(0, 0, 1, 12) - AIC:2779.4242140348188
Current function value: 22.737751
Iterations: 2
Function evaluations: 45
Gradient evaluations: 34
SARIMA(1, 1, 3)x(0, 0, 3, 12) - AIC:7377.123261323514
Optimization terminated successfully.
Current function value: 6.027724
Iterations: 176
Function evaluations: 380
Gradient evaluations: 380
SARIMA(1, 1, 2)x(2, 0, 1, 12) - AIC:2473.016953321524
Optimization terminated successfully.
Current function value: 6.686578
Iterations: 120
Function evaluations: 237
Gradient evaluations: 237
SARIMA(1, 1, 3)x(1, 0, 0, 12) - AIC:2727.3134137372526
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 6.292448
Iterations: 88
Function evaluations: 184
Gradient evaluations: 184
SARIMA(1, 1, 3)x(0, 0, 2, 12) - AIC:2433.76419472674
SARIMA(1, 1, 3)x(1, 0, 3, 12) - AIC:7692.477574980119
Current function value: 5.943157
Iterations: 200
Function evaluations: 422
Gradient evaluations: 422
SARIMA(1, 1, 2)x(1, 0, 2, 12) - AIC:2439.46792951787
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 3)x(2, 0, 3, 12) - AIC:7724.96927538047
Current function value: 5.471164
Iterations: 200
Function evaluations: 442
Gradient evaluations: 442
SARIMA(0, 1, 3)x(3, 0, 2, 12) - AIC:2244.5552183054942
Current function value: 5.934711
Iterations: 200
Function evaluations: 428
Gradient evaluations: 428
SARIMA(1, 1, 2)x(2, 0, 2, 12) - AIC:2431.5016362434335
Current function value: 6.453730
Iterations: 200
Function evaluations: 423
Gradient evaluations: 422
SARIMA(1, 1, 3)x(1, 0, 1, 12) - AIC:2644.5549983502765
Optimization terminated successfully.
Current function value: 5.433299
Iterations: 109
Function evaluations: 239
Gradient evaluations: 239
SARIMA(1, 1, 2)x(3, 0, 0, 12) - AIC:2230.605452582383
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 3)x(3, 0, 3, 12) - AIC:7788.132354495019
Optimization terminated successfully.
Current function value: 9.408458
Iterations: 31
Function evaluations: 94
Gradient evaluations: 94
SARIMA(2, 1, 0)x(0, 0, 0, 12) - AIC:3867.598458606961
Optimization terminated successfully.
Current function value: 7.717429
Iterations: 29
Function evaluations: 51
Gradient evaluations: 51
SARIMA(2, 1, 0)x(0, 0, 1, 12) - AIC:3026.349089524774
Optimization terminated successfully.
Current function value: 6.036413
Iterations: 166
Function evaluations: 330
Gradient evaluations: 330
SARIMA(1, 1, 3)x(2, 0, 0, 12) - AIC:2471.382922022125
Current function value: 25.929190
Iterations: 0
Function evaluations: 18
Gradient evaluations: 6
SARIMA(2, 1, 0)x(0, 0, 3, 12) - AIC:8238.21016494751
Optimization terminated successfully.
Current function value: 7.067944
Iterations: 58
Function evaluations: 126
Gradient evaluations: 126
SARIMA(2, 1, 0)x(1, 0, 0, 12) - AIC:3151.265090611391
Optimization terminated successfully.
Current function value: 6.685069
Iterations: 35
Function evaluations: 75
Gradient evaluations: 75
SARIMA(2, 1, 0)x(0, 0, 2, 12) - AIC:2599.891602938661
Optimization terminated successfully.
Current function value: 6.736934
Iterations: 38
Function evaluations: 73
Gradient evaluations: 73
SARIMA(2, 1, 0)x(1, 0, 1, 12) - AIC:2740.4389152449394
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 0)x(1, 0, 3, 12) - AIC:7866.645190810615
Optimization terminated successfully.
Current function value: 6.125165
Iterations: 113
Function evaluations: 230
Gradient evaluations: 230
SARIMA(2, 1, 0)x(1, 0, 2, 12) - AIC:2499.7549247862703
Optimization terminated successfully.
Current function value: 6.090282
Iterations: 91
Function evaluations: 165
Gradient evaluations: 165
SARIMA(2, 1, 0)x(2, 0, 1, 12) - AIC:2485.079845455093
Current function value: 20.566248
Iterations: 200
Function evaluations: 406
Gradient evaluations: 406
SARIMA(2, 1, 0)x(2, 0, 0, 12) - AIC:2208.075713813471
Current function value: 6.019251
Iterations: 200
Function evaluations: 431
Gradient evaluations: 431
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(1, 1, 3)x(2, 0, 1, 12) - AIC:2453.0905375370608
SARIMA(2, 1, 0)x(2, 0, 3, 12) - AIC:7900.024998989915
Current function value: 5.862805
Iterations: 200
Function evaluations: 413
Gradient evaluations: 413
SARIMA(1, 1, 3)x(1, 0, 2, 12) - AIC:2379.8370538858276
Current function value: 5.431710
Iterations: 200
Function evaluations: 412
Gradient evaluations: 412
SARIMA(1, 1, 1)x(3, 0, 2, 12) - AIC:2230.6530285195377
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 0)x(3, 0, 3, 12) - AIC:7880.342897381581
Optimization terminated successfully.
Current function value: 8.390111
Iterations: 32
Function evaluations: 63
Gradient evaluations: 63
SARIMA(2, 1, 1)x(0, 0, 0, 12) - AIC:3290.350547642971
Optimization terminated successfully.
Current function value: 7.313005
Iterations: 47
Function evaluations: 107
Gradient evaluations: 107
SARIMA(2, 1, 1)x(0, 0, 1, 12) - AIC:2853.2017227820256
Optimization terminated successfully.
Current function value: 6.483665
Iterations: 40
Function evaluations: 79
Gradient evaluations: 79
Optimization terminated successfully.
Current function value: 5.431842
Iterations: 170
Function evaluations: 356
Gradient evaluations: 356
SARIMA(2, 1, 1)x(0, 0, 2, 12) - AIC:2499.6500694223246
SARIMA(1, 1, 2)x(3, 0, 1, 12) - AIC:2231.969501449708
Current function value: 5.858315
Iterations: 200
Function evaluations: 429
Gradient evaluations: 429
SARIMA(1, 1, 3)x(2, 0, 2, 12) - AIC:2403.9808571136155
Optimization terminated successfully.
Current function value: 6.750113
Iterations: 39
Function evaluations: 70
Gradient evaluations: 70
SARIMA(2, 1, 1)x(1, 0, 0, 12) - AIC:2745.444357154154
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 1)x(1, 0, 3, 12) - AIC:7800.710722289244
Optimization terminated successfully.
Current function value: 6.650956
Iterations: 61
Function evaluations: 135
Gradient evaluations: 135
SARIMA(2, 1, 1)x(1, 0, 1, 12) - AIC:2684.4863797142043
Current function value: 23.228059
Iterations: 1
Function evaluations: 73
Gradient evaluations: 61
SARIMA(2, 1, 1)x(0, 0, 3, 12) - AIC:8163.3428203781505
Optimization terminated successfully.
Current function value: 6.090756
Iterations: 39
Function evaluations: 76
Gradient evaluations: 76
SARIMA(2, 1, 1)x(2, 0, 0, 12) - AIC:2484.0705489295106
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 1)x(2, 0, 3, 12) - AIC:7833.794494541582
Optimization terminated successfully.
Current function value: 5.465385
Iterations: 73
Function evaluations: 155
Gradient evaluations: 155
SARIMA(2, 1, 0)x(3, 0, 1, 12) - AIC:2233.7948656579556
Current function value: 6.076659
Iterations: 200
Function evaluations: 419
Gradient evaluations: 419
SARIMA(2, 1, 0)x(2, 0, 2, 12) - AIC:2482.0047830266
Optimization terminated successfully.
Current function value: 5.580438
Iterations: 191
Function evaluations: 384
Gradient evaluations: 384
SARIMA(2, 1, 0)x(3, 0, 0, 12) - AIC:2611.201373202886
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 1)x(3, 0, 3, 12) - AIC:7898.145206260745
Current function value: 5.433260
Iterations: 200
Function evaluations: 433
Gradient evaluations: 433
SARIMA(1, 1, 2)x(3, 0, 2, 12) - AIC:2235.2778065060643
Optimization terminated successfully.
Current function value: 8.223054
Iterations: 116
Function evaluations: 238
Gradient evaluations: 238
SARIMA(2, 1, 2)x(0, 0, 0, 12) - AIC:3238.012013475489
Optimization terminated successfully.
Current function value: 5.462515
Iterations: 83
Function evaluations: 162
Gradient evaluations: 162
SARIMA(2, 1, 1)x(3, 0, 0, 12) - AIC:2231.573767186974
Current function value: 25.467475
Iterations: 0
Function evaluations: 9
Gradient evaluations: 6
SARIMA(2, 1, 2)x(0, 0, 3, 12) - AIC:8087.590222744408
Optimization terminated successfully.
Current function value: 5.972775
Iterations: 167
Function evaluations: 349
Gradient evaluations: 349
SARIMA(2, 1, 1)x(1, 0, 2, 12) - AIC:2450.1015166759967
Optimization terminated successfully.
Current function value: 7.168999
Iterations: 72
Function evaluations: 143
Gradient evaluations: 143
SARIMA(2, 1, 2)x(0, 0, 1, 12) - AIC:2806.095468004935
Current function value: 5.421539
Iterations: 200
Function evaluations: 426
Gradient evaluations: 426
SARIMA(1, 1, 3)x(3, 0, 0, 12) - AIC:2223.5622141013655
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 2)x(1, 0, 3, 12) - AIC:7686.275125516444
Optimization terminated successfully.
Current function value: 6.630137
Iterations: 101
Function evaluations: 189
Gradient evaluations: 189
SARIMA(2, 1, 2)x(1, 0, 0, 12) - AIC:2714.1889075838226
Optimization terminated successfully.
Current function value: 6.356640
Iterations: 78
Function evaluations: 170
Gradient evaluations: 170
SARIMA(2, 1, 2)x(0, 0, 2, 12) - AIC:2463.542309987522
Current function value: 6.060077
Iterations: 200
Function evaluations: 436
Gradient evaluations: 436
SARIMA(2, 1, 1)x(2, 0, 1, 12) - AIC:2468.166277272711
Optimization terminated successfully.
Current function value: 6.500506
Iterations: 93
Function evaluations: 214
Gradient evaluations: 214
SARIMA(2, 1, 2)x(1, 0, 1, 12) - AIC:2644.2471796082336
Optimization terminated successfully.
Current function value: 6.053355
Iterations: 177
Function evaluations: 365
Gradient evaluations: 365
SARIMA(2, 1, 1)x(2, 0, 2, 12) - AIC:2474.7964918142566
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 2)x(2, 0, 3, 12) - AIC:7719.06286184443
Current function value: 5.419112
Iterations: 200
Function evaluations: 433
Gradient evaluations: 433
SARIMA(1, 1, 3)x(3, 0, 1, 12) - AIC:2223.5508453775024
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 2)x(3, 0, 3, 12) - AIC:7782.8197572596455
Optimization terminated successfully.
Current function value: 5.990571
Iterations: 166
Function evaluations: 332
Gradient evaluations: 332
SARIMA(2, 1, 2)x(2, 0, 0, 12) - AIC:2458.6003102702102
Current function value: 8.124487
Iterations: 200
Function evaluations: 449
Gradient evaluations: 449
SARIMA(2, 1, 3)x(0, 0, 0, 12) - AIC:3231.9280879463163
Current function value: 5.416913
Iterations: 200
Function evaluations: 451
Gradient evaluations: 451
SARIMA(1, 1, 3)x(3, 0, 2, 12) - AIC:2218.81176350927
Current function value: 5.464554
Iterations: 200
Function evaluations: 429
Gradient evaluations: 429
SARIMA(2, 1, 0)x(3, 0, 2, 12) - AIC:2235.1526612944613
Current function value: 22.730281
Iterations: 0
Function evaluations: 29
Gradient evaluations: 19
SARIMA(2, 1, 3)x(0, 0, 3, 12) - AIC:7999.451622551792
Current function value: 5.924858
Iterations: 200
Function evaluations: 382
Gradient evaluations: 382
SARIMA(2, 1, 2)x(1, 0, 2, 12) - AIC:2435.5238860059962
Optimization terminated successfully.
Current function value: 5.458925
Iterations: 170
Function evaluations: 346
Gradient evaluations: 346
SARIMA(2, 1, 1)x(3, 0, 1, 12) - AIC:2230.7179565103806
Optimization terminated successfully.
Current function value: 5.967994
Iterations: 150
Function evaluations: 348
Gradient evaluations: 348
SARIMA(2, 1, 2)x(2, 0, 1, 12) - AIC:2433.8859778664964
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 3)x(1, 0, 3, 12) - AIC:7698.34311441394
Optimization terminated successfully.
Current function value: 7.109604
Iterations: 176
Function evaluations: 362
Gradient evaluations: 362
SARIMA(2, 1, 3)x(0, 0, 1, 12) - AIC:2780.9977107541113
Optimization terminated successfully.
Current function value: 6.625064
Iterations: 191
Function evaluations: 354
Gradient evaluations: 354
SARIMA(2, 1, 3)x(1, 0, 0, 12) - AIC:2707.8175902149287
Current function value: 5.909883
Iterations: 200
Function evaluations: 419
Gradient evaluations: 419
SARIMA(2, 1, 2)x(2, 0, 2, 12) - AIC:2385.1022960228333
Optimization terminated successfully.
Current function value: 6.291338
Iterations: 116
Function evaluations: 245
Gradient evaluations: 245
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(2, 1, 3)x(0, 0, 2, 12) - AIC:2435.7769393430817
SARIMA(2, 1, 3)x(2, 0, 3, 12) - AIC:7730.834814814844
Current function value: 6.451884
Iterations: 200
Function evaluations: 434
Gradient evaluations: 434
SARIMA(2, 1, 3)x(1, 0, 1, 12) - AIC:2620.8158994989562
Current function value: 5.372639
Iterations: 200
Function evaluations: 403
Gradient evaluations: 403
SARIMA(2, 1, 1)x(3, 0, 2, 12) - AIC:2208.778034023946
Optimization terminated successfully.
Current function value: 8.354967
Iterations: 82
Function evaluations: 184
Gradient evaluations: 184
SARIMA(3, 1, 0)x(0, 0, 0, 12) - AIC:4162.750408496155
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 7.704669
Iterations: 17
Function evaluations: 30
Gradient evaluations: 30
SARIMA(3, 1, 0)x(0, 0, 1, 12) - AIC:3022.5396920609264
SARIMA(2, 1, 3)x(3, 0, 3, 12) - AIC:7793.99789392884
Optimization terminated successfully.
Current function value: 5.375358
Iterations: 170
Function evaluations: 353
Gradient evaluations: 353
SARIMA(2, 1, 2)x(3, 0, 0, 12) - AIC:2209.562191080884
Current function value: 25.915756
Iterations: 1
Function evaluations: 25
Gradient evaluations: 16
SARIMA(3, 1, 0)x(0, 0, 3, 12) - AIC:8237.148429181409
Optimization terminated successfully.
Current function value: 6.677827
Iterations: 56
Function evaluations: 102
Gradient evaluations: 102
SARIMA(3, 1, 0)x(0, 0, 2, 12) - AIC:2597.879082584913
Current function value: 5.987355
Iterations: 200
Function evaluations: 393
Gradient evaluations: 393
SARIMA(2, 1, 3)x(2, 0, 0, 12) - AIC:2452.722368454258
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 0)x(1, 0, 3, 12) - AIC:7889.805572597701
Optimization terminated successfully.
Current function value: 6.747779
Iterations: 198
Function evaluations: 407
Gradient evaluations: 407
SARIMA(3, 1, 0)x(1, 0, 0, 12) - AIC:3227.010625505511
Optimization terminated successfully.
Current function value: 6.653125
Iterations: 126
Function evaluations: 290
Gradient evaluations: 290
SARIMA(3, 1, 0)x(1, 0, 1, 12) - AIC:2710.7003978697207
Optimization terminated successfully.
Current function value: 6.106297
Iterations: 95
Function evaluations: 192
Gradient evaluations: 192
SARIMA(3, 1, 0)x(1, 0, 2, 12) - AIC:2494.7466765676527
Current function value: 5.970674
Iterations: 200
Function evaluations: 428
Gradient evaluations: 428
SARIMA(2, 1, 3)x(2, 0, 1, 12) - AIC:2449.5628737875804
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 0)x(2, 0, 3, 12) - AIC:7923.185380777448
Optimization terminated successfully.
Current function value: 6.024587
Iterations: 82
Function evaluations: 164
Gradient evaluations: 164
SARIMA(3, 1, 0)x(2, 0, 1, 12) - AIC:2460.7000535416505
Current function value: 5.862398
Iterations: 200
Function evaluations: 407
Gradient evaluations: 407
SARIMA(2, 1, 3)x(1, 0, 2, 12) - AIC:2381.4568952392547
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Current function value: 12.270641
Iterations: 200
Function evaluations: 420
Gradient evaluations: 420
SARIMA(3, 1, 0)x(2, 0, 0, 12) - AIC:2354.23492276185
SARIMA(3, 1, 0)x(3, 0, 3, 12) - AIC:7819.91429413181
Optimization terminated successfully.
Current function value: 8.319739
Iterations: 29
Function evaluations: 72
Gradient evaluations: 72
SARIMA(3, 1, 1)x(0, 0, 0, 12) - AIC:3266.1142638141155
Optimization terminated successfully.
Current function value: 7.305688
Iterations: 64
Function evaluations: 121
Gradient evaluations: 121
SARIMA(3, 1, 1)x(0, 0, 1, 12) - AIC:2852.9568690184133
Current function value: 5.372280
Iterations: 200
Function evaluations: 420
Gradient evaluations: 420
SARIMA(2, 1, 2)x(3, 0, 1, 12) - AIC:2207.1967791985903
Current function value: 25.682718
Iterations: 0
Function evaluations: 17
Gradient evaluations: 7
SARIMA(3, 1, 1)x(0, 0, 3, 12) - AIC:8166.545331313209
Optimization terminated successfully.
Current function value: 6.669261
Iterations: 36
Function evaluations: 66
Gradient evaluations: 66
SARIMA(3, 1, 1)x(1, 0, 0, 12) - AIC:2718.4348216854182
Current function value: 5.857852
Iterations: 200
Function evaluations: 408
Gradient evaluations: 408
SARIMA(2, 1, 3)x(2, 0, 2, 12) - AIC:2387.533454673543
Optimization terminated successfully.
Current function value: 6.468627
Iterations: 77
Function evaluations: 156
Gradient evaluations: 156
SARIMA(3, 1, 1)x(0, 0, 2, 12) - AIC:2498.5024780390904
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 1)x(1, 0, 3, 12) - AIC:7802.710722289244
Current function value: 5.372560
Iterations: 200
Function evaluations: 417
Gradient evaluations: 417
SARIMA(2, 1, 2)x(3, 0, 2, 12) - AIC:2213.1280365327216
Optimization terminated successfully.
Current function value: 6.578271
Iterations: 126
Function evaluations: 251
Gradient evaluations: 251
SARIMA(3, 1, 1)x(1, 0, 1, 12) - AIC:2660.554759348346
Optimization terminated successfully.
Current function value: 6.024378
Iterations: 77
Function evaluations: 145
Gradient evaluations: 145
SARIMA(3, 1, 1)x(2, 0, 0, 12) - AIC:2460.874006565255
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Current function value: 6.007612
Iterations: 200
Function evaluations: 415
Gradient evaluations: 415
SARIMA(3, 1, 0)x(2, 0, 2, 12) - AIC:2454.2361458285686
SARIMA(3, 1, 1)x(2, 0, 3, 12) - AIC:7835.794494541976
Current function value: 5.372557
Iterations: 200
Function evaluations: 387
Gradient evaluations: 387
SARIMA(2, 1, 3)x(3, 0, 0, 12) - AIC:2204.9378079401595
Optimization terminated successfully.
Current function value: 5.972265
Iterations: 159
Function evaluations: 322
Gradient evaluations: 322
SARIMA(3, 1, 1)x(1, 0, 2, 12) - AIC:2451.944220601348
Current function value: 22.789291
Iterations: 200
Function evaluations: 422
Gradient evaluations: 422
SARIMA(3, 1, 0)x(3, 0, 0, 12) - AIC:1898.2904638910143
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 1)x(3, 0, 3, 12) - AIC:7815.750243641939
Optimization terminated successfully.
Current function value: 5.394266
Iterations: 118
Function evaluations: 233
Gradient evaluations: 233
SARIMA(3, 1, 0)x(3, 0, 1, 12) - AIC:2207.807494653924
Optimization terminated successfully.
Current function value: 8.221896
Iterations: 99
Function evaluations: 203
Gradient evaluations: 203
SARIMA(3, 1, 2)x(0, 0, 0, 12) - AIC:3238.7223263551105
Optimization terminated successfully.
Current function value: 5.936193
Iterations: 166
Function evaluations: 337
Gradient evaluations: 337
SARIMA(3, 1, 1)x(2, 0, 1, 12) - AIC:2436.0567871394105
Current function value: 24.231569
Iterations: 0
Function evaluations: 32
Gradient evaluations: 27
SARIMA(3, 1, 2)x(0, 0, 3, 12) - AIC:7157.345250168024
Optimization terminated successfully.
Current function value: 7.166134
Iterations: 111
Function evaluations: 200
Gradient evaluations: 200
SARIMA(3, 1, 2)x(0, 0, 1, 12) - AIC:2806.662175187309
Optimization terminated successfully.
Current function value: 5.392334
Iterations: 96
Function evaluations: 195
Gradient evaluations: 195
SARIMA(3, 1, 1)x(3, 0, 0, 12) - AIC:2207.272210946737
Optimization terminated successfully.
Current function value: 6.578005
Iterations: 87
Function evaluations: 159
Gradient evaluations: 159
SARIMA(3, 1, 2)x(1, 0, 0, 12) - AIC:2693.993486604333
Optimization terminated successfully.
Current function value: 6.350516
Iterations: 59
Function evaluations: 121
Gradient evaluations: 121
SARIMA(3, 1, 2)x(0, 0, 2, 12) - AIC:2462.0382754444386
Current function value: 25.126930
Iterations: 0
Function evaluations: 23
Gradient evaluations: 11
SARIMA(3, 1, 2)x(1, 0, 3, 12) - AIC:8000.353723396716
Current function value: 25.233832
Iterations: 0
Function evaluations: 24
Gradient evaluations: 13
SARIMA(3, 1, 2)x(2, 0, 3, 12) - AIC:8033.141355686761
Optimization terminated successfully.
Current function value: 6.500783
Iterations: 76
Function evaluations: 181
Gradient evaluations: 181
SARIMA(3, 1, 2)x(1, 0, 1, 12) - AIC:2672.1380627030744
Current function value: 5.370284
Iterations: 200
Function evaluations: 403
Gradient evaluations: 403
SARIMA(2, 1, 3)x(3, 0, 1, 12) - AIC:2205.674346270177
Current function value: 5.923987
Iterations: 200
Function evaluations: 409
Gradient evaluations: 409
SARIMA(3, 1, 1)x(2, 0, 2, 12) - AIC:2433.1997285409375
Optimization terminated successfully.
Current function value: 5.941913
Iterations: 93
Function evaluations: 189
Gradient evaluations: 189
SARIMA(3, 1, 2)x(2, 0, 0, 12) - AIC:2439.6252414367586
Current function value: 25.448267
Iterations: 0
Function evaluations: 40
Gradient evaluations: 29
SARIMA(3, 1, 2)x(3, 0, 3, 12) - AIC:8096.898281806342
Optimization terminated successfully.
Current function value: 8.166604
Iterations: 122
Function evaluations: 278
Gradient evaluations: 278
SARIMA(3, 1, 3)x(0, 0, 0, 12) - AIC:3207.8662693890733
Optimization terminated successfully.
Current function value: 5.910366
Iterations: 128
Function evaluations: 291
Gradient evaluations: 291
SARIMA(3, 1, 2)x(2, 0, 2, 12) - AIC:2411.113693677575
Current function value: 22.476365
Iterations: 1
Function evaluations: 23
Gradient evaluations: 13
SARIMA(3, 1, 3)x(0, 0, 3, 12) - AIC:7688.416006224009
Current function value: 5.369332
Iterations: 200
Function evaluations: 451
Gradient evaluations: 451
Optimization terminated successfully.
Current function value: 5.917288
Iterations: 197
Function evaluations: 398
Gradient evaluations: 398
SARIMA(2, 1, 3)x(3, 0, 2, 12) - AIC:2193.0387069369435
SARIMA(3, 1, 2)x(1, 0, 2, 12) - AIC:2406.8729360479665
Current function value: 5.392798
Iterations: 200
Function evaluations: 418
Gradient evaluations: 418
SARIMA(3, 1, 0)x(3, 0, 2, 12) - AIC:2209.465068609356
Current function value: 7.090432
Iterations: 200
Function evaluations: 436
Gradient evaluations: 436
SARIMA(3, 1, 3)x(0, 0, 1, 12) - AIC:2771.8469781670556
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 3)x(1, 0, 3, 12) - AIC:7739.9075414487725
Current function value: 5.921582
Iterations: 200
Function evaluations: 432
Gradient evaluations: 432
SARIMA(3, 1, 2)x(2, 0, 1, 12) - AIC:2416.7150160538686
Current function value: 6.569849
Iterations: 200
Function evaluations: 337
Gradient evaluations: 337
SARIMA(3, 1, 3)x(1, 0, 0, 12) - AIC:2683.281283493975
Optimization terminated successfully.
Current function value: 6.280406
Iterations: 160
Function evaluations: 290
Gradient evaluations: 290
SARIMA(3, 1, 3)x(0, 0, 2, 12) - AIC:2512.3264121565394
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
SARIMA(3, 1, 3)x(2, 0, 3, 12) - AIC:7772.402148920781
Current function value: 6.450198
Iterations: 200
Function evaluations: 422
Gradient evaluations: 422
SARIMA(3, 1, 3)x(1, 0, 1, 12) - AIC:2618.496690345624
Optimization terminated successfully.
Current function value: 5.318886
Iterations: 145
Function evaluations: 310
Gradient evaluations: 310
SARIMA(3, 1, 2)x(3, 0, 0, 12) - AIC:2188.4399312834994
Current function value: 5.317594
Iterations: 200
Function evaluations: 438
Gradient evaluations: 438
SARIMA(3, 1, 1)x(3, 0, 1, 12) - AIC:2184.456427882312
Optimization terminated successfully.
Current function value: -0.000000
Iterations: 1
Function evaluations: 13
Gradient evaluations: 13
Optimization terminated successfully.
Current function value: 5.937401
Iterations: 151
Function evaluations: 286
Gradient evaluations: 286
SARIMA(3, 1, 3)x(3, 0, 3, 12) - AIC:7835.562241335231
SARIMA(3, 1, 3)x(2, 0, 0, 12) - AIC:2430.03763499995
Current function value: 5.316546
Iterations: 200
Function evaluations: 409
Gradient evaluations: 409
SARIMA(3, 1, 1)x(3, 0, 2, 12) - AIC:2187.1213601035524
Current function value: 5.859664
Iterations: 200
Function evaluations: 407
Gradient evaluations: 407
SARIMA(3, 1, 3)x(1, 0, 2, 12) - AIC:2379.6322681119923
Current function value: 5.915791
Iterations: 200
Function evaluations: 429
Gradient evaluations: 429
SARIMA(3, 1, 3)x(2, 0, 1, 12) - AIC:2411.03783547713
Current function value: 5.855423
Iterations: 200
Function evaluations: 394
Gradient evaluations: 394
SARIMA(3, 1, 3)x(2, 0, 2, 12) - AIC:2385.522317710321
Optimization terminated successfully.
Current function value: 5.315724
Iterations: 199
Function evaluations: 424
Gradient evaluations: 424
SARIMA(3, 1, 2)x(3, 0, 1, 12) - AIC:2188.905101053825
Current function value: 5.312900
Iterations: 200
Function evaluations: 434
Gradient evaluations: 434
SARIMA(3, 1, 2)x(3, 0, 2, 12) - AIC:2183.2145796593927
Current function value: 5.315229
Iterations: 200
Function evaluations: 387
Gradient evaluations: 387
SARIMA(3, 1, 3)x(3, 0, 0, 12) - AIC:2180.901169851145
Current function value: 5.312806
Iterations: 200
Function evaluations: 416
Gradient evaluations: 416
SARIMA(3, 1, 3)x(3, 0, 1, 12) - AIC:2182.0342780061933
Current function value: 5.310372
Iterations: 200
Function evaluations: 420
Gradient evaluations: 420
SARIMA(3, 1, 3)x(3, 0, 2, 12) - AIC:2182.4513420867515
SARIMA_AIC.sort_values(by='AIC',ascending=True)
| param | seasonal | AIC | |
|---|---|---|---|
| 213 | (3, 1, 0) | (3, 0, 0, 12) | 1898.290464 |
| 253 | (3, 1, 3) | (3, 0, 0, 12) | 2180.901170 |
| 254 | (3, 1, 3) | (3, 0, 1, 12) | 2182.034278 |
| 255 | (3, 1, 3) | (3, 0, 2, 12) | 2182.451342 |
| 252 | (3, 1, 2) | (3, 0, 2, 12) | 2183.214580 |
| ... | ... | ... | ... |
| 73 | (1, 1, 1) | (0, 0, 3, 12) | 8167.707947 |
| 2 | (0, 1, 0) | (0, 0, 3, 12) | 8235.438068 |
| 185 | (3, 1, 0) | (0, 0, 3, 12) | 8237.148429 |
| 121 | (2, 1, 0) | (0, 0, 3, 12) | 8238.210165 |
| 60 | (1, 1, 0) | (0, 0, 3, 12) | 8240.146956 |
256 rows × 3 columns
###fir the model
import statsmodels.api as sm
##i have had to use cfg processing and it takes quite a wile
auto_SARIMA = sm.tsa.statespace.SARIMAX(train_data['Sparkling_Clipped'],order=(3,1,0),seasonal_order=(3,0,0,12),enforce_stationarity=False,
enforce_invertibility=False)
results_auto_SARIMA = auto_SARIMA.fit(method = "cg")
print(results_auto_SARIMA.summary())
Current function value: 33.315046
Iterations: 50
Function evaluations: 105
Gradient evaluations: 105
SARIMAX Results
==========================================================================================
Dep. Variable: Sparkling_Clipped No. Observations: 144
Model: SARIMAX(3, 1, 0)x(3, 0, 0, 12) Log Likelihood -900.204
Date: Sun, 29 Oct 2023 AIC 1814.407
Time: 19:29:05 BIC 1832.918
Sample: 01-31-1980 HQIC 1821.907
- 12-31-1991
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.4538 2.284 -0.199 0.843 -4.931 4.023
ar.L2 -0.5227 3.009 -0.174 0.862 -6.421 5.375
ar.L3 -0.1308 3.229 -0.041 0.968 -6.460 6.199
ar.S.L12 0.5404 2.570 0.210 0.833 -4.496 5.577
ar.S.L24 0.2222 3.009 0.074 0.941 -5.675 6.120
ar.S.L36 0.2670 2.799 0.095 0.924 -5.218 5.752
sigma2 5.06e+06 1.05e+06 4.824 0.000 3e+06 7.12e+06
===================================================================================
Ljung-Box (L1) (Q): 1.50 Jarque-Bera (JB): 25.37
Prob(Q): 0.22 Prob(JB): 0.00
Heteroskedasticity (H): 0.80 Skew: 0.52
Prob(H) (two-sided): 0.51 Kurtosis: 5.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
###running hte diagnostics
results_auto_SARIMA.plot_diagnostics();
There are a few outliers on the Q-Q plot but everything else looks in line with the theoretical expectations and there is no pattern to be derived from the plots
#Predict on the Test Set using this model and evaluate the model.
predicted_auto_SARIMA = results_auto_SARIMA.get_forecast(steps=len(test_data))
#get the head of the data
predicted_auto_SARIMA.summary_frame(alpha=0.05).head()
| Sparkling_Clipped | mean | mean_se | mean_ci_lower | mean_ci_upper |
|---|---|---|---|---|
| 1992-01-31 | 1717.166653 | 2249.546691 | -2691.863843 | 6126.197149 |
| 1992-02-29 | 1259.098096 | 2563.218262 | -3764.717383 | 6282.913574 |
| 1992-03-31 | 1476.780926 | 2614.680619 | -3647.898919 | 6601.460771 |
| 1992-04-30 | 1179.226095 | 2828.514581 | -4364.560614 | 6723.012804 |
| 1992-05-31 | 1213.039611 | 3125.360808 | -4912.555011 | 7338.634234 |
##Get the rmse and caluclate the error
rmse = mean_squared_error(test_data['Sparkling'],(predicted_auto_SARIMA.predicted_mean),squared=False)
mape = mean_absolute_percentage_error(test_data['Sparkling'],(predicted_auto_SARIMA.predicted_mean))
print('RMSE:',rmse,'\nMAPE:',mape)
RMSE: 704.493874007528 MAPE: 29.56426350152451
##add resukts to the data frame
resultsDf_SARIMA = pd.DataFrame({'Test RMSE': [sqrt(mean_squared_error(test_data['Sparkling'],(predicted_auto_SARIMA.predicted_mean)))]}
,index=['SARIMA(3,1,0)(3,0,0,12)'])
resultsDf = pd.concat([resultsDf, resultsDf_SARIMA])
resultsDf
| Test RMSE | |
|---|---|
| RegressionOnTime | 1344.644564 |
| NaiveModel | 3979.914692 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
| Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing | 304.152583 |
| ARIMA(2,1,2) | 1300.857871 |
| SARIMA(0,1,1)(1,0,1,12) | 1225.638727 |
| SARIMA(3,1,0)(3,0,0,12) | 704.493874 |
temp_resultsDf = pd.DataFrame({'RMSE': rmse,'MAPE':mape}
,index=['SARIMA(3,1,0)(3,0,0,12)'])
plt.plot(train_data['Sparkling'], label='Train')
plt.plot(test_data['Sparkling'], label='Test')
plt.plot((predicted_auto_SARIMA.predicted_mean), label='SARIMA(3,1,0)(3,0,0,12) predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("SARIMA Forecast")
plt.legend(loc='best')
plt.grid();
In this section we choose which model performs the best on the test data set
###sort the results
resultsDf.sort_values(by = 'Test RMSE')
| Test RMSE | |
|---|---|
| Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing | 304.152583 |
| auto_TES, Alpha=0.077016,Beta=0.077016,Gamma=0.306 | 351.507911 |
| SARIMA(3,1,0)(3,0,0,12) | 704.493874 |
| SimpleAverageModel | 1268.462755 |
| auto_SES,Alpha =0.0370 | 1290.978200 |
| ARIMA(2,1,2) | 1300.857871 |
| RegressionOnTime | 1344.644564 |
| Alpha=0.1,SimpleExponentialSmoothing | 1356.185563 |
| Alpha=0.1,Beta=0.1,DoubleExponentialSmoothing | 1357.769783 |
| auto_DES, Alpha = 0.737,Beta = 0.0001 | 1992.838366 |
| NaiveModel | 3979.914692 |
plt.plot(df_final['Sparkling'], label='Actual Dataset')
plt.plot(TES_test['predict',index_best_model], label='Alpha=0.1,Beta=0.2,Gamma=0.2,TripleExponentialSmoothing predictions on Test Set')
plt.plot(TES_test['auto_predict'], label='Alpha=0.077016,Beta=0.077016,Gamma=0.306,Auto TripleExponentialSmoothing predictions on Test Set')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("Top two models Forecast")
plt.legend(loc='best')
plt.grid();
The TES models are the best performing models
The TES model found through hyper-parameter selection is the best perfroming model
The TES model has the following parameters:
We have selected the most optimum model for the data and we will proceed to build the model and predict 12 months into the future.
####see the results from the TES Model
df_results_TES.sort_values(by =['Test RMSE'], ascending = True).head()
| index | trend | seasonal_periods | seasonal | smoothing_level | smoothing_trend | smoothing_seasonal | damped_trend | damping_trend | Train RMSE | Test RMSE | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 17517 | 17517 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.4 | 385.434751 | 304.152583 |
| 16059 | 16059 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.3 | 385.434751 | 304.152583 |
| 13143 | 13143 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.1 | 385.434751 | 304.152583 |
| 24807 | 24807 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.9 | 385.434751 | 304.152583 |
| 18975 | 18975 | mul | 12 | mul | 0.1 | 0.2 | 0.2 | False | 0.5 | 385.434751 | 304.152583 |
###getting the best model variables from the results table
trend = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][1]
seasonal_periods = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][2]
seasonal = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][3]
smoothing_level = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][4]
smoothing_trend = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][5]
smoothing_seasonal = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][6]
damped_trend = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][7]
damping_trend = df_results_TES.sort_values(by =['Test RMSE'], ascending = True).values[0][8]
####Setting the final model
final_model= ExponentialSmoothing(TES_train['Sparkling_Clipped'], trend=trend,
damped_trend=damped_trend,
seasonal_periods=seasonal_periods, freq = 'M',
seasonal=seasonal
#,use_boxcox=use_boxcox
).fit(
# use_boxcox = use_boxcox,
smoothing_level=smoothing_level,
smoothing_trend=smoothing_trend,
smoothing_seasonal= smoothing_seasonal,
damping_trend=damping_trend,
optimized=False,
use_brute=True)
##setting the predictions for the year
predicitons_final = final_model.predict(start='31-08-1995',end='31-08-1996')
plt.plot(df_final['Sparkling'], label='Actual Dataset')
plt.plot(predicitons_final,label='Alpha =0.1, Beta = Gamma = 0.2, Triple Exponential Smoothing Model')
plt.xlabel('Year')
plt.ylabel('Sparkling')
plt.title("TES Final Forecast")
plt.legend(loc='best')
plt.grid();
####calculating the prediction interval
pred_1_df = pd.DataFrame({'lower_CI':predicitons_final - 1.96*np.std(final_model.resid,ddof=1),
'prediction':predicitons_final,
'upper_ci':predicitons_final + 1.96*np.std(final_model.resid,ddof=1)})
pred_1_df.head()
| lower_CI | prediction | upper_ci | |
|---|---|---|---|
| 1995-08-31 | 1466.439234 | 2160.616141 | 2854.793049 |
| 1995-09-30 | 1709.045668 | 2403.222575 | 3097.399483 |
| 1995-10-31 | 2603.338646 | 3297.515553 | 3991.692461 |
| 1995-11-30 | 3633.749013 | 4327.925920 | 5022.102827 |
| 1995-12-31 | 5346.591176 | 6040.768083 | 6734.944990 |
##look at the tail of the data
pred_1_df.tail()
| lower_CI | prediction | upper_ci | |
|---|---|---|---|
| 1996-04-30 | 1038.892794 | 1733.069701 | 2427.246609 |
| 1996-05-31 | 942.991674 | 1637.168581 | 2331.345488 |
| 1996-06-30 | 821.131776 | 1515.308684 | 2209.485591 |
| 1996-07-31 | 1404.631833 | 2098.808741 | 2792.985648 |
| 1996-08-31 | 1485.479788 | 2179.656695 | 2873.833603 |
# plot the forecast along with the confidence band
axis = df_final["Sparkling"].plot(label='Actual')
pred_1_df['prediction'].plot(ax=axis, label='Forecast', alpha=0.5)
axis.fill_between(pred_1_df.index, pred_1_df['lower_CI'], pred_1_df['upper_ci'], color='k', alpha=.15)
axis.set_xlabel('Year-Months')
axis.set_ylabel('Sparkling')
plt.title("Final Forecast")
plt.legend(loc='best')
plt.grid()
plt.show()
We will look at a few things:
#Sales in the first years
print(f""" We have sold {sum(df_final['Sparkling'])} of Sparkling to date""")
We have sold 449252 of Sparkling to date
###the sales for the previous 12 months
print(f""" In the last 12 months we sold {sum(df_final["Sparkling"].tail(12))} of Sparkling""")
In the last 12 months we sold 29196 of Sparkling
## sales in the next 12 month
print(f""" We are expected to sell {sum(predicitons_final)} of Sparkling inb the next 12 month""")
We are expected to sell 32765.92244576343 of Sparkling inb the next 12 month
### look at the how the predicted sales compare to the last 12 months
YOY_CHANGE = ((sum(predicitons_final) - sum(df_final["Sparkling"].tail(12)) ) / sum(df_final["Sparkling"].tail(12)) ) * 100
print(f"""We will sell {YOY_CHANGE}% more Sparkling in the next 12 months compared to the previous 12 months""")
We will sell 12.227436791901049% more Sparkling in the next 12 months compared to the previous 12 months
quarterly_prediction = predicitons_final.groupby(predicitons_final.index.quarter).sum().values
quarterly_prediction
array([ 5371.86177046, 4885.54696624, 8842.30415281, 13666.20955626])
quarterly_actuals= df_final["Sparkling"].tail(12).groupby(df_final["Sparkling"].tail(12).index.quarter).sum().values
quarterly_actuals
array([ 4369, 5220, 6494, 13113])
##percentage change quarter on quarter
((quarterly_prediction - quarterly_actuals)/quarterly_actuals) * 100
array([22.95403457, -6.40714624, 36.16113571, 4.21878713])
The growth in the sales is driven by stronger expected sales for the next 12 months in quarter 1 and quarter 3 compared to the previous 12 months. We expect to see a 36.2% growth in sales for quarter 3 and a 23% growth in sales for quarter one.We also see a 4% growth in quarter four sales. The business is expected to grow quarter on quarter except for quarter 2.